{ "cells": [ { "metadata": {}, "cell_type": "markdown", "source": "# Chapter 6: Markov Chain Monte Carlo", "id": "56f5bbe4a34bf86a" }, { "metadata": {}, "cell_type": "markdown", "source": [ "## Introduction to Markov chains\n", "Markov chain Monte Carlo (Markov Chain simulation or MCMC for short) is a Bayesian way of finding values of $\\theta$ from an approximate distribution, where $\\theta$ are the unknown parameters that make up the distribution. The Markov chain algorithm corrects this approximate distribution in each iteration, eventually converging to a final distribution that aligns with the posterior distribution $p(\\theta|y)$ (for more on the posterior distribution, see [Chapter 4](https://bayesian-statistics-for-astrophysics-2024.readthedocs.io/en/latest/lecture_notes/group4/Chapter4.html)). This is particularly useful when the posterior is difficult or impossible to sample making MCMC a powerful tool. This converging of the distribution is where the power of Markov chain simulations lies, allowing for the determination of the original distribution solely from its data. \n", "\n", "This is done by constructing a stationary Markov chain, a sequence of random variables, $\\theta_t$, where $\\theta_t$ has the same distribution as the posterior. Now instead of taking independent variables, the Markov chain introduces variables $\\theta_{\\{1, ..., T\\}}$ obeying the Markov property\n", "\n", "$$\n", "\tp(\\theta_{t+1} | \\theta_{t},...,\\theta_1) = p(\\theta_{t+1}|\\theta_{t}),\n", "$$\n", "such that the variable $\\theta_{t+1}$ is only dependent on the most recent variable $\\theta_{t}$ and therefore, the distribution of $\\theta_{t+1}$ does not depend on how the chain got to $\\theta_{t}$, but only what this state is. In other words, one could say that \"_What comes next is determined entirely by the present situation._\"\n", "\n", "Markov chain simulations are frequently used when it is inefficient to find $\\theta$ directly from the posterior $p(\\theta|y).$ As mentioned before, we sample a new, random $\\theta$ iteratively, with the expectation that the distribution converges more and more to $p(\\theta|y).$ Generally, MCMC is exceptionally effective for high-dimensional distributions, as direct-sampling Monte Carlo methods cannot complete such analysis in a reasonable time. \n", "\n", "MCMC algorithms, such as Metropolis-Hasting and affine parameters, propose and accept new steps of the MCMC chain in the parameter space, allowing the chain to explore the parameter space." ], "id": "61675086c38cd5d5" }, { "metadata": {}, "cell_type": "markdown", "source": [ "## Explanation of Metropolis-Hasting\n", "The simplest method of simulating a Markov chain is the Metropolis-Hastings algorithm, which itself relies on the more basic Metropolis algorithm. This we will inquire about first.\n", "\n", "The Metropolis algorithm expands on a random walk, a process where each step is determined randomly based on the current position, with a rule that accepts or rejects the sampled data to converge the sampled data to the specified target distribution. This accepting and rejecting is based on a given rule. The algorithm works as follows:\n", "\n", "1. Draw a starting point $\\theta_0$ with $p(\\theta_0|y) > 0$ from a certain starting distribution $p_0(\\theta)$. For example, $p_0(\\theta) \\sim \\mathcal{N}(\\theta_{t-1}, \\sigma^2)$ if $\\theta \\in \\mathbb{R}$.\n", "\n", "2. For $t = 1, 2, \\ldots$:\n", "\n", " a. Sample a proposal $\\theta_\\ast$ from a proposal distribution $J_t(\\theta_\\ast|\\theta_{t-1})$ at time $t$. We require $J_t(\\theta_a|\\theta_b) = J_t(\\theta_b|\\theta_a)$ for all $\\theta_a$, $\\theta_b$, and $t$ for now.\n", "\n", " b. Calculate the ratio of the densities of $\\theta_a$ and $\\theta_b$:\n", "\n", " $$\n", " r = \\frac{p(\\theta_a|y)}{p(\\theta_b|y)} = \\frac{p(\\theta_{\\ast}|y)}{p(\\theta_{t-1}|y)}.\n", " $$\n", " \n", " c. Set \n", "\n", " $$\n", " \\theta_t = \\begin{cases}\n", " \\theta_\\ast &\\ \\text{if } u \\text{ } \\leq \\min(r, 1), \\text{where } u \\sim \\text{Uniform(0,1)} \\\\\n", " \\theta_{t-1} &\\ \\text{otherwise}.\n", " \\end{cases}\n", " $$\n", "\n", "To put this into words; a value of $\\theta_*$ is sampled from a previous distribution, and if the ratio of densities $r$ is lower than a random uniform variable, the proposed state is rejected and $\\theta_t$ = $\\theta_{t-1}$.\n", "\n", "Now we return to the Metropolis-Hastings (MH) algorithm. In contrast to the Metropolis algorithm, the MH algorithm does not require $J_t(\\theta_*|\\theta_{t-1})$ to be symmetric. Due to this new, asymmetric proposal distribution, we must alter the ratio of densities of $\\theta_a$ and $\\theta_b$, which is as follows:\n", "\n", "$$ \n", "r = \\frac{p(\\theta_{\\ast}|y)J_t(\\theta_{t-1}|\\theta_{\\ast})}{p(\\theta_{t-1}| y)J_t(\\theta_{\\ast}|\\theta_{t-1})}\n", "$$\n", "\n", "For the proposal distribution, the following characteristics are favorable. First, it is easy to sample from $J(\\theta_\\ast| \\theta)$ for any $\\theta$. Second, the ratio of $r$ (as given by the above equation) is easy to calculate. Third, the random walk does not go too slowly, i.e., the sampled values are not too close together. Finally, rejected jumps do not happen often, as they hurt the performance of the algorithm.\n", "\n", "The acceptance rate of the MH algorithm is the proportion of the steps that are accepted. A good MH algorithm has in general an acceptance rate of 0.44. When the acceptance rate is too high (e.g, >70\\%), the MH algorithm is likely taking very small steps resulting in a slow exploration of parameter space. When the acceptance rate is too low (e.g, <10\\%), the MH algorithm is rejecting too many steps, making it inefficient.\n" ], "id": "773203b39c8c2ff2" }, { "metadata": {}, "cell_type": "markdown", "source": [ "### Burn-in Time and Autocorrelation\n", "In the Metropolis-Hastings algorithm, the chain takes some steps to reach a stable value, called the burn-in time. We discard these points before we make the probability distribution. In a one-parameter fit, this can easily be done by eye. However, in fits with many parameters, you would prefer an automated process to take care of this. This can be done via autocorrelation. If we have the chain $\\{X_1, X_2,...,X_n\\}$ we can calculate the autocorrelation coefficient with\n", "$$\n", "\\rho(k)=\\frac{1}{(n-k)\\sigma^2}\\sum_{t=1}^{n-k}(X_t-\\mu)(X_{t+k}-\\mu) \\tag{6.1}\n", "$$\n", "where $k$ is the lag or time gap between the samples and $\\mu$ the mean. We can choose a threshold value, say 0.1, and calculate above which value of $k$ the autocorrelation coefficient is below the threshold. We can then take $\\{X_k,...,X_n\\}$ to be samples from our posterior distribution. In the example below we will use autocorrelation to determine the burn-in time." ], "id": "10efc5f499bc238c" }, { "metadata": {}, "cell_type": "markdown", "source": [ "### 1D Example\n", "We have collected data and want to infer the mean and its error. The data points are $y_i=$[1,2,3] with error $\\sigma_{y_i}=$[1,1,1]. We can do this analytically and find $\\mu=2$ and $\\sigma_\\mu=1/\\sqrt{3}\\approx0.577$. We want to derive these values using the Metropolis-Hasting algorithm. First, we need to construct the posterior. We will assume a flat prior making the posterior the likelihood. Next, we pick a starting point for the Markov chain. Using this as the center of a normal distribution, a new point is drawn from the distribution. We find the ratio of the new points posterior to the last point. If the ratio is larger than a random value between 0 and 1, we accept the new point to the chain. Otherwise, the last point value is added to the chain. These steps are repeated until we can reconstruct the desired distribution for $\\mu$. We take 10 as the starting point and use a starting distribution with a scale parameter of 3. This allows the chain to converge to the correct value quickly. \n", "\n", "As mentioned before, the chain takes some steps to get to a stable value. We disregard these points before we make the probability distribution. \n", "\n" ], "id": "9d3743ce2d883a33" }, { "metadata": { "ExecuteTime": { "end_time": "2024-12-19T18:26:33.280260Z", "start_time": "2024-12-19T18:26:28.452679Z" } }, "cell_type": "code", "source": [ "import numpy as np\n", "import matplotlib.pyplot as plt\n", "from scipy.stats import norm\n", "\n", "steps = 10000\n", "X = 10\n", "Xs = []\n", "y = np.array([1, 2, 3])\n", "y_sigma = np.array([1, 1, 1])\n", "\n", "\n", "def prob_dist(mu):\n", " return np.exp(-1 * np.sum(((y - mu) ** 2 / (2 * y_sigma ** 2))))\n", "\n", "\n", "def autocorrelation(chain, lag):\n", " sum = np.sum((chain[:len(chain) - lag] - np.mean(chain)) * (chain[lag:] - np.mean(chain)))\n", " normalization = (len(chain) - lag) * np.var(chain)\n", " autocorr = sum / normalization\n", " return autocorr\n", "\n", "\n", "for _ in range(steps):\n", " Y = np.random.normal(X, 1)\n", " q = prob_dist(Y) / prob_dist(X)\n", " if q >= np.random.rand():\n", " X = Y\n", " Xs.append(X)\n", "\n", "auto = np.inf\n", "index = 1\n", "while auto >= 0.01:\n", " auto = autocorrelation(Xs, index)\n", " index += 1\n", " if index >= len(Xs):\n", " print('autocorrelation failed')\n", " break\n", "\n", "print(f'The burn-in time is at step {index}.')\n", "plt.figure()\n", "plt.plot(Xs[:10 * index])\n", "plt.vlines(index, 0, 10, colors='red', linestyles='dashed')\n", "plt.xlabel('Step')\n", "plt.ylabel(r'$\\mu$')\n", "plt.title('Simulated MCMC with burn-in time.')\n", "plt.show()\n", "\n", "plt.figure()\n", "plt.hist(Xs[index:], bins=steps // 100, density=True)\n", "mu = np.mean(Xs[index:])\n", "mu_sigma = np.std(Xs[index:])\n", "print(f'Simulated mean: {mu:.3f}, True mean: 2')\n", "print(f'Simulated error: {mu_sigma:.3f}, True error: 0.577')\n", "x = np.linspace(mu - 4 * mu_sigma, mu + 4 * mu_sigma, 1000)\n", "plt.plot(x, norm.pdf(x, mu, mu_sigma))\n", "plt.xlabel(r'$\\mu$')\n", "plt.ylabel('Probability Density')\n", "plt.title(r'Simulated distribution of $\\mu$ after burn-in.')\n", "plt.show()" ], "id": "464a4e5c214c98a4", "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "The burn-in time is at step 15.\n" ] }, { "data": { "text/plain": [ "
" ], "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjMAAAHHCAYAAABKudlQAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy80BEi2AAAACXBIWXMAAA9hAAAPYQGoP6dpAABrs0lEQVR4nO3dd3xT9f4/8NdJmiad6Z50sWRvRIaigqLXgQuvXlRArxcVB/ITt6LXwZXrQL2K617Ar6CIVxS9igOZsqfMstpS6KKlbTqTJjm/P07PadOmbdqmTU76ej4eeZhmfo5kvPP+vD/vjyCKoggiIiIildJ4egBERERE7cFghoiIiFSNwQwRERGpGoMZIiIiUjUGM0RERKRqDGaIiIhI1RjMEBERkaoxmCEiIiJVYzBDREREqsZghlQjNTUV06dP98hzv/DCCxAEwSPP3ZTMzEwIgoAlS5Z4eiiqNX36dKSmprp82+Dg4DY/V2pqKq699to237+zXXrppbj00ks75bk8+d4m38BghjzuwIEDuOWWW5CSkgKDwYDExERcccUVePfddz09NLd4//33PRpwrF+/HoIgQBAEfPbZZ05vM3bsWAiCgAEDBjS6zmazYfHixbj00ksREREBvV6P1NRUzJgxA7t27VJut2TJEuV5Nm/e3OhxRFFEUlISBEFw+qVeXV2Nt956C6NGjYLRaITBYEDv3r3x4IMP4tixY+34P+C6yspKvPDCC1i/fn2nPF9XsmXLFrzwwgsoKSnx9FDIB/l5egDUtW3ZsgWXXXYZkpOTce+99yIuLg7Z2dnYtm0b3n77bTz00EPKbdPT06HRqC/+fv/99xEVFeXxX54GgwHLly/HHXfc4XB5ZmYmtmzZAoPB0Og+VVVVuOmmm7BmzRpccsklePrppxEREYHMzEx8+eWXWLp0KU6fPo1u3bo1ep5x48Y5PNaGDRtw5swZ6PX6Rs9TWFiIq666Crt378a1116Lv/zlLwgODkZ6ejq++OILfPTRR7BYLG76P1Hn448/ht1uV/6urKzEiy++CACdlpXwVj///LNbH2/Lli148cUXMX36dISFhTlcp9b3NnkPBjPkUa+88gqMRiN27tzZ6AOuoKDA4W9nX4Lkuj/96U9YvXo1CgsLERUVpVy+fPlyxMbGolevXiguLna4z9y5c7FmzRq89dZbmD17tsN18+bNw1tvveX0eVauXIl33nkHfn51HzHLly/H8OHDUVhY2Og+06dPx969e/HVV1/h5ptvdrjupZdewjPPPNOWQ26RTqfrkMf1FFEUUV1djYCAgHY/lr+/vxtG5Bq+t6m9GAqTR508eRL9+/dvFMgAQExMjMPfDefV5WmNzZs34+GHH0Z0dDTCwsIwc+ZMWCwWlJSU4K677kJ4eDjCw8Px+OOPo/4m8fL0S8MpBVdrURYvXozLL78cMTEx0Ov16NevHxYtWtRozIcOHcKGDRuUKZj6v/hLSkowe/ZsJCUlQa/Xo2fPnnjttdccsgXy7aZPnw6j0YiwsDBMmzat1en6yZMnQ6/XY+XKlQ6XL1++HLfeeiu0Wq3D5WfOnMGHH36IK664olEgAwBarRaPPfaYQ1YGAG6//XYUFRXhl19+US6zWCz46quv8Je//KXR42zfvh3/+9//cM899zQKZADpi+71119v8rhKSkqg1WrxzjvvKJcVFhZCo9EgMjLS4d/8/vvvR1xcnPJ3/ZqZzMxMREdHAwBefPFF5d/rhRdecHi+s2fP4oYbbkBwcDCio6Px2GOPwWazNTm+hn7++WcMGTIEBoMB/fr1w9dff+1wfVP1WfLrPTMzU7lMrsP56aefMGLECAQEBODDDz9UXttffvklXnnlFXTr1g0GgwETJkzAiRMnXBpnw5qZ9jzmCy+8gLlz5wIA0tLSlP+38rG4+70NAHa7HQsXLkT//v1hMBgQGxuLmTNnNgrYyTcwM0MelZKSgq1bt+LgwYNO6zVc8dBDDyEuLg4vvvgitm3bho8++ghhYWHYsmULkpOT8eqrr+KHH37AP//5TwwYMAB33XWXW8a+aNEi9O/fH9dffz38/Pzw3Xff4YEHHoDdbsesWbMAAAsXLsRDDz2E4OBgJbsQGxsLQJrSGD9+PM6ePYuZM2ciOTkZW7ZswVNPPYXc3FwsXLgQgPRre/Lkydi8eTPuu+8+9O3bF6tWrcK0adNaNd7AwEBMnjwZn3/+Oe6//34AwP79+3Ho0CF88skn+OOPPxxu/+OPP8JqteLOO+9s1fOkpqZi9OjR+Pzzz3H11Vcrj1VaWorbbrvNIegAgNWrVwNAq59HFhYWhgEDBmDjxo14+OGHAQCbN2+GIAg4f/48Dh8+jP79+wMANm3ahIsvvtjp40RHR2PRokW4//77ceONN+Kmm24CAAwaNEi5jc1mw6RJkzBq1Ci8/vrr+PXXX/HGG2+gR48eyv/T5hw/fhx//vOfcd9992HatGlYvHgxpkyZgjVr1uCKK65o0/Gnp6fj9ttvx8yZM3HvvffiggsuUK77xz/+AY1Gg8ceewylpaVYsGABpk6diu3bt7fpudr6mDfddBOOHTuGzz//HG+99ZaSGZSDx6a05709c+ZMLFmyBDNmzMDDDz+MjIwM/Otf/8LevXvx+++/+1xWrssTiTzo559/FrVarajVasXRo0eLjz/+uPjTTz+JFoul0W1TUlLEadOmKX8vXrxYBCBOmjRJtNvtyuWjR48WBUEQ77vvPuUyq9UqduvWTRw/frxy2bp160QA4rp16xyeJyMjQwQgLl68WLls3rx5YsO3S2VlZaMxTpo0SezevbvDZf3793d4XtlLL70kBgUFiceOHXO4/MknnxS1Wq14+vRpURRF8ZtvvhEBiAsWLHA4nosvvrjROJ2Rj3PlypXi999/LwqCoDz23LlzlfGOHz9e7N+/v3K/Rx99VAQg7t27t9nHl8n/Hjt37hT/9a9/iSEhIcr/oylTpoiXXXaZKIrSv+M111yj3O/GG28UAYjFxcUuPY8zs2bNEmNjY5W/58yZI15yySViTEyMuGjRIlEURbGoqEgUBEF8++23ldtNmzZNTElJUf4+d+6cCECcN29eo+eYNm2aCED8+9//7nD50KFDxeHDh7c4xpSUFBGA+N///le5rLS0VIyPjxeHDh2qXObstSaKdf9/MzIyGj3mmjVrHG4r/5v37dtXNJvNyuVvv/22CEA8cOBAi+MdP3680/dLWx/zn//8Z6Px1z8Od763N23aJAIQly1b5vA8a9ascXo5qR+nmcijrrjiCmzduhXXX3899u/fjwULFmDSpElITExUfrG35J577nFIy48aNQqiKOKee+5RLtNqtRgxYgROnTrltrHXr0soLS1FYWEhxo8fj1OnTqG0tLTF+69cuRIXX3wxwsPDUVhYqJwmTpwIm82GjRs3AgB++OEH+Pn5Ofzy12q1DsXRrrryyisRERGBL774AqIo4osvvsDtt9/u9LYmkwkAEBIS0urnufXWW1FVVYXvv/8eZWVl+P77751OMbX3eWQXX3wx8vPzkZ6eDkDKwFxyySW4+OKLsWnTJgBStkYUxSYzM6667777Gj23q6+rhIQE3HjjjcrfoaGhuOuuu7B3717k5eW1aTxpaWmYNGmS0+tmzJjhUPsiH3t73gcd8ZhNaet7e+XKlTAajbjiiisc3lvDhw9HcHAw1q1b5/axkmdxmok8buTIkfj6669hsViwf/9+rFq1Cm+99RZuueUW7Nu3D/369Wv2/snJyQ5/G41GAEBSUlKjy905X/77779j3rx52Lp1KyorKx2uKy0tVcbRlOPHj+OPP/5oMtUuF0BnZWUhPj6+UY+T+tMJrtLpdJgyZQqWL1+OCy+8ENnZ2U0GGaGhoQCAsrKyVj9PdHQ0Jk6ciOXLl6OyshI2mw233HJLi8/jrHbKFfIX6qZNm9CtWzfs3bsXL7/8MqKjo5V6m02bNiE0NBSDBw9u03MA0kqthv9e4eHhLr+uevbs2agepnfv3gCkmp369TyuSktLa/K6hu+N8PBwAFDGW15ejvLycuV6rVbb4tRPS4/pTm19bx8/fhylpaWN6u5kDRcXkPoxmCGv4e/vj5EjR2LkyJHo3bs3ZsyYgZUrV2LevHnN3q9h4Wpzl4v1igSbaoLnSjHnyZMnMWHCBPTp0wdvvvkmkpKS4O/vjx9++AFvvfVWowJeZ+x2O6644go8/vjjTq+Xv+Tc7S9/+Qs++OADvPDCCxg8eHCTwWKfPn0ASH2AhgwZ0qbnuffee5GXl4err766yUCl/vO0NWuSkJCAtLQ0bNy4EampqRBFEaNHj0Z0dDQeeeQRZGVlYdOmTRgzZky7lgA39Vpzp9a+LptbudTUeOX3weuvv64sRQekGrb6BcZteUx3aut72263IyYmBsuWLXN6/5YCNlIfBjPklUaMGAEAyM3N7bDnkH9RNlwVlJWV1eJ9v/vuO5jNZqxevdrh16Oz9HVTX049evRAeXk5Jk6c2OxzpaSkYO3atSgvL3fIzshTKq01btw4JCcnY/369XjttdeavN3VV18NrVaLzz77rE3FuTfeeCNmzpyJbdu2YcWKFU3e7rrrrsP8+fPx2WeftWsK6OKLL8bGjRuRlpaGIUOGICQkBIMHD4bRaMSaNWuwZ88ehy9uZzq6y/OJEycgiqLD88gNAeVVVfVfl/UDQFdel6111113OfQDcseS7qZ0ZgftHj164Ndff8XYsWM79JjIe7Bmhjxq3bp1Tn/R/fDDDwDaNpXiqpSUFGi1WqU2Rfb++++3eF/5l2H9sZeWlmLx4sWNbhsUFOR0GfWtt96KrVu34qeffmp0XUlJCaxWKwCpb4vVanVY9m2z2drcIVkQBLzzzjuYN29es0FKUlIS7r33Xvz8889On8tut+ONN97AmTNnnN4/ODgYixYtwgsvvIDrrruuyecZPXo0rrrqKnzyySf45ptvGl1vsVjw2GOPtXhcF198MTIzM7FixQolKNJoNBgzZgzefPNN1NTUtBgsBQYGAmgc4LpLTk4OVq1apfxtMpnw6aefYsiQIcoUU48ePQDA4XVZUVGBpUuXun083bt3x8SJE5XT2LFj3f4csqCgIAAd9/+2vltvvRU2mw0vvfRSo+usVqvDGHJzc3H06FHU1NR0+Lio4zAzQx710EMPobKyEjfeeCP69OkDi8WCLVu2YMWKFUrL/I5iNBoxZcoUvPvuuxAEAT169MD333/v0nz6lVdeCX9/f1x33XWYOXMmysvL8fHHHyMmJqZRNmn48OFYtGgRXn75ZfTs2RMxMTG4/PLLMXfuXKxevRrXXnstpk+fjuHDh6OiogIHDhzAV199hczMTERFReG6667D2LFj8eSTTyIzM1PpTeJKkXFTJk+ejMmTJ7d4uzfeeAMnT57Eww8/jK+//hrXXnstwsPDcfr0aaxcuRJHjx7Fbbfd1uT9XV0+/umnn+LKK6/ETTfdhOuuuw4TJkxAUFAQjh8/ji+++AK5ubnN9poB6upm0tPT8eqrryqXX3LJJfjxxx+h1+sxcuTIZh8jICAA/fr1w4oVK9C7d29ERERgwIABbW4b0FDv3r1xzz33YOfOnYiNjcV//vMf5OfnOwTBV155JZKTk3HPPfdg7ty50Gq1+M9//oPo6GicPn3aLePwhOHDhwMAnnnmGdx2223Q6XS47rrrlCDHncaPH4+ZM2di/vz52LdvH6688krodDocP34cK1euxNtvv63UcD311FNYunQpMjIyXN6ni7wPgxnyqNdffx0rV67EDz/8oLSsT05OxgMPPIBnn322zQWhrnr33XdRU1ODDz74AHq9HrfeeqvSs6I5F1xwAb766is8++yzeOyxxxAXF4f7778f0dHRuPvuux1u+/zzzyMrKwsLFixAWVkZxo8fj8svvxyBgYHYsGEDXn31VaxcuRKffvopQkND0bt3b7z44otKsaNGo8Hq1asxe/ZsfPbZZxAEAddffz3eeOMNDB06tMP+3wBSpuLHH3/EkiVLsHTpUrz00kuorKxEQkICLr/8cixbtgyJiYntfp7o6Ghs2bIF77//PlasWIFnnnkGFosFKSkpuP766/HII4+0+BgXXHABYmJiUFBQ4DB1Igc5F154oUudZj/55BM89NBDePTRR2GxWDBv3jy3BTO9evXCu+++i7lz5yI9PR1paWlYsWKFw2oknU6HVatW4YEHHsBzzz2HuLg4zJ49G+Hh4R0a3He0kSNH4qWXXsIHH3yANWvWwG63IyMjo0OCGQD44IMPMHz4cHz44Yd4+umn4efnh9TUVNxxxx0dmoEizxDEjqjaIiIiIuokrJkhIiIiVWMwQ0RERKrGYIaIiIhUjcEMERERqRqDGSIiIlI1BjNERESkal2iz4zdbkdOTg5CQkI6taU2ERERtZ0oiigrK0NCQkKz+6p1iWAmJyen0S6rREREpA7Z2dno1q1bk9d3iWAmJCQEgPQ/IzQ01MOjISIiIleYTCYkJSUp3+NN6RLBjDy1FBoaymCGiIhIZVoqEWEBMBEREakagxkiIiJSNQYzREREpGoMZoiIiEjVGMwQERGRqjGYISIiIlVjMENERESqxmCGiIiIVI3BDBEREakagxkiIiJSNY8HMxs3bsR1112HhIQECIKAb775xuF6URTx/PPPIz4+HgEBAZg4cSKOHz/umcESERGR1/F4MFNRUYHBgwfjvffec3r9ggUL8M477+CDDz7A9u3bERQUhEmTJqG6urqTR0pERETeyOMbTV599dW4+uqrnV4niiIWLlyIZ599FpMnTwYAfPrpp4iNjcU333yD2267rTOHSkRERF7I45mZ5mRkZCAvLw8TJ05ULjMajRg1ahS2bt3a5P3MZjNMJpPDiYiIiHyTVwczeXl5AIDY2FiHy2NjY5XrnJk/fz6MRqNySkpK6tBxqkpFBSAI0qmiwtOjISIiajevDmba6qmnnkJpaalyys7O9vSQiIiIqIN4dTATFxcHAMjPz3e4PD8/X7nOGb1ej9DQUIcTERER+SavDmbS0tIQFxeHtWvXKpeZTCZs374do0eP9uDIiIiIyFt4fDVTeXk5Tpw4ofydkZGBffv2ISIiAsnJyZg9ezZefvll9OrVC2lpaXjuueeQkJCAG264wXODVjONBhg/vu48ERGRynk8mNm1axcuu+wy5e85c+YAAKZNm4YlS5bg8ccfR0VFBf72t7+hpKQE48aNw5o1a2AwGDw1ZHULCADWr/f0KIiIiNxGEEVR9PQgOprJZILRaERpaSnrZ4iIiFTC1e9vzjMQERGRqjGY6WoqKoDoaOnEPjNEROQDPF4zQx5QWOjpERAREbkNMzNERESkagxmiIiISNUYzBAREZGqMZghIiIiVWMwQ0RERKrG1UxdjUYDjBhRd56IiEjlGMx0NQEBwM6dnh4FERGR2/CnOREREakagxkiIiJSNQYzXU1lJZCaKp0qKz09GiIionZjzUxXI4pAVlbdeSIiIpVjZoaIiIhUjcEMERERqRqDGSIiIlI1BjNERESkagxmiIiISNW4mqmrEQSgX7+680RERCrHYKYd1h7Jx4GzpcrfwXo/TBmRBGOAzoOjakFgIHDokKdHQURE5DYMZtph7dECLN9+2uGy6hobHry8l4dGRERE1PUwmGmHUWkR0NTO1BzKMWHv6RKcPs+uukRERJ2JwUw7TB6SiMlDEgEAK3aext7TJSgoM3t4VC2orARGjpTO79wpTTsRERGpGIMZN4kJMQAACkxeHsyIInD4cN15IiIilePSbDeJDtEDAM6Ve3kwQ0RE5GMYzLhJTG0wU1Ruhs3OjAcREVFnYTDjJpHBemgEwC5KAQ0RERF1DgYzbqLVCIgMlrIzXl8ETERE5EMYzLiRPNVUUFbt4ZEQERF1HVzN5EZKEbA3Z2YEAUhJqTtPRESkcgxm3EjJzHjz8uzAQCAz09OjICIichtOM7mR0mvGmzMzREREPobBjBvFhLJmhoiIqLMxmHGjGDXUzFRVSdsZjBwpnSciIlI51sy4UbQappnsdmDXrrrzREREKsfMjBvVLc02Q+S+R0RERJ2CwYwbyUuzLVY7TFVWD4+GiIioa2Aw40YGnRahBmnmjkXAREREnYPBjJvFhEp1M15dBExERORDGMy4Wf26GSIiIup4XM3kZtFq2J8pKsrTIyAiInIbBjNu5vVbGgQFAefOeXoUREREbsNpJjeTtzQ4V+6lwQwREZGPYTDjZsqWBt6amSEiIvIxDGbczOtrZqqqgEsvlU7czoCIiHwAa2bczOtXM9ntwIYNdeeJiIhUjpkZN5P3ZyqrtqK6xubh0RAREfk+BjNuFmrwg95P+t/KxnlEREQdj8GMmwmCUFcE7K11M0RERD6EwUwHkJdnc0UTERFRx2MBcAeIDpYyMzszixEW6A+dVsCARCMMOq2HR0ZEROR7GMx0gNjaaab//J6B//yeAUCqpbl2cAJuGJKoXK/TahBvNEAQhM4dYGBg5z4fERFRB2Iw0wFuHt4N+7JLUGGRVjOVVFpQWG7B8u2nsXz7aYfbPnZlbzx4ea/OG1xQEFBR0XnPR0RE1MEEURRFTw+io5lMJhiNRpSWliI0NLTTn99mF7HtVBH+u+cM1qefg7nGBrsIVNXYEBXsj61PTYBOy/IlIiKi+lz9/mZmphNoNQLG9ozC2J51u1XX2OwYPf83FJabsfZIAa4aEOfBERIREamX16cDbDYbnnvuOaSlpSEgIAA9evTASy+9BLUnlHRaDW4enggA+HJXduc9cXU1cM010qmaS8eJiEj9vD4z89prr2HRokVYunQp+vfvj127dmHGjBkwGo14+OGHPT28dvnziCR8uOEU1qcXIK+0GnFGQ8c/qc0G/PBD3XkiIiKV8/rMzJYtWzB58mRcc801SE1NxS233IIrr7wSO3bs8PTQ2q17dDAuTI2AXQS+2t2J2RkiIiIf4vXBzJgxY7B27VocO3YMALB//35s3rwZV199dZP3MZvNMJlMDidv9eeRSQCAFbuyYbere+qMiIjIE7w+mHnyySdx2223oU+fPtDpdBg6dChmz56NqVOnNnmf+fPnw2g0KqekpKROHHHr/GlgPEL0fsg+X4Vtp4o8PRwiIiLV8fpg5ssvv8SyZcuwfPly7NmzB0uXLsXrr7+OpUuXNnmfp556CqWlpcopO9t7p3AC/LW4bkgCAGD1/hwPj4aIiEh9vL4AeO7cuUp2BgAGDhyIrKwszJ8/H9OmTXN6H71eD71e35nDbJcLUyOwfPtpZBaxmR0REVFreX1mprKyEhqN4zC1Wi3sdruHRuR+8bWrmHJLuVSaiIiotbw+M3PdddfhlVdeQXJyMvr374+9e/fizTffxN133+3poblNQlgAACmYEUWxY/dqCgoCVN6jh4iIqD6vD2beffddPPfcc3jggQdQUFCAhIQEzJw5E88//7ynh+Y2saEGCAJgsdpRVGFBVLB6psiIiIg8zeuDmZCQECxcuBALFy709FA6jL+fBtHBehSUmZFbUs1ghoiIqBW8vmamq4ivnWrKKa3q2CeqrgamTJFO3M6AiIh8AIMZL5FQWwScU9LBwYzNBnz1lXTidgZEROQDGMx4iXhjXREwERERuY7BjJdICOukzAwREZGPYTDjJZiZISIiahsGM14ivjYzk8vMDBERUaswmPESibWrmfLLzLBx92wiIiKXMZjxElHBevhpBNjsIgrKONVERETkKgYzXkKrERAb2glFwIGBQHm5dAoM7LjnISIi6iQMZrxI3YqmDszMCIK0P1NQkHSeiIhI5RjMeJG6FU0sAiYiInIVgxkvEt8ZmRmzGZg+XTqZzR33PERERJ2EwYwXSeiMzIzVCixdKp2s1o57HiIiok7CYMaLJISxcR4REVFrMZjxIvHGTphmIiIi8jEMZryInJkpLDfDbOWO1kRERK5gMONFwgN10PtJ/yR5nGoiIiJyCYMZLyIIgpKd4VQTERGRaxjMeBm5boa9ZoiIiFzj5+kBkKO6xnkdlJkJDAQKCurOExERqRyDGS+TWNs477NtWdh47Bx0Wg0euKwHxvSIcs8TCAIQHe2exyIiIvICDGa8TN/4UABSZkbOzmg1gvuCGSIiIh/DYMbLTOofh+V/HYXiyhoczTPh3d9OoKDMjdsOmM3AnDnS+TffBPR69z02ERGRB7AA2MtoNALG9IzCNYPiMal/HACp74zbWK3A++9LJ25nQEREPoDBjBeLDpGyJkXlZtjsoodHQ0RE5J0YzHixiCB/CAJgF4HiSounh0NEROSVGMx4MZ1Wg/BAfwBunmoiIiLyIQxmvFxUcG0wU8bMDBERkTMMZrxcVLBUN8PMDBERkXMMZrycXAR8zp3Ls4mIiHwI+8x4ObdnZgICgIyMuvNEREQqx2DGy8nBzDl3BTMaDZCa6p7HIiIi8gKcZvJySgFwOQuAiYiInGEw4+XcXjNjsQBz50onCwMkIiJSPwYzXs7tNTM1NcDrr0unmhr3PCYREZEHMZjxcnJm5nyFBXZuaUBERNQIgxkvFxEk1czY7CK3NCAiInKCwYyXk7Y00AFw44omIiIiH8JgRgXkqSZuaUBERNQYgxkV4JYGRERETWMwowIMZoiIiJrGDsAq4NYuwAEBwMGDdeeJiIhUjsGMCri1cZ5GA/Tv3/7HISIi8hKcZlIBbmlARETUNGZmVCBKWc3khsyMxQK8+qp0/umnAX//9j8mERGRBzGYUYFodxYA19QAL74onZ87l8EMERGpHqeZVEAuAC7ilgZERESNMJhRgchgbmlARETUFAYzKlB/SwMWARMRETliMKMSbJxHRETkHIMZlVAa57ljRRMREZEPYTCjEspmk8zMEBEROeDSbJVw25YGBgOwY0fdeSIiIpVjMKMSUSHSiqY9WcVYsfM0tBoNLr0gWglyXKbVAiNHdsAIiYiIPIPBjErEhUpZlJ2ZxdiZWQwAuKp/HD64c7gnh0VERORxDGZU4op+sbhtZBIKy80oqrBg7+kSZBdXtv6BLBbg7bel8488wg7ARESkeqooAD579izuuOMOREZGIiAgAAMHDsSuXbs8PaxOFWLQ4R83D8In00bimT/1BQCUVVtb/0A1NcDjj0unmho3j5KIiKjzeX1mpri4GGPHjsVll12GH3/8EdHR0Th+/DjCw8M9PTSPCTFIDfTKqhmMEBEReX0w89prryEpKQmLFy9WLktLS/PgiDwvxCD9s5VVWyGKIgRB8PCIiIiIPMfrp5lWr16NESNGYMqUKYiJicHQoUPx8ccfN3sfs9kMk8nkcPIlcjBjtYuorrF7eDRERESe5fXBzKlTp7Bo0SL06tULP/30E+6//348/PDDWLp0aZP3mT9/PoxGo3JKSkrqxBF3vCB/P8jJGE41ERFRVyeIoih6ehDN8ff3x4gRI7Blyxblsocffhg7d+7E1q1bnd7HbDbDbK5rLmcymZCUlITS0lKEhoZ2+Jg7w8AXfkJZtRW/zhmPnjHBrt+xogIIrr19eTkQFNQxAyQiImonk8kEo9HY4ve312dm4uPj0a9fP4fL+vbti9OnTzd5H71ej9DQUIeTrwllETAREREAFRQAjx07Funp6Q6XHTt2DCkpKR4akXeoXwTcKgYDsG5d3XkiIiKV8/pg5tFHH8WYMWPw6quv4tZbb8WOHTvw0Ucf4aOPPvL00DyqzcGMVgtceqn7B0REROQhXj/NNHLkSKxatQqff/45BgwYgJdeegkLFy7E1KlTPT00j2KvGSIiIonXZ2YA4Nprr8W1117r6WF4lTZnZmpqADmr9be/ATqdm0dGRETUuVQRzFBjdcFMKzMzFgvw4IPS+enTGcwQEZHqef00EzknTzOZ2rI/ExERkQ9hMKNSbZ5mIiIi8jEMZlSKBcBEREQSBjMqFcrMDBEREQAGM6qlTDOZmZkhIqKujcGMStVNMzEzQ0REXRuXZqtUmwuA9Xrg++/rzhMREakcgxmVql8ALIoiBEFw7Y5+fsA113TgyIiIiDoXp5lUSs7M1NhEmK12D4+GiIjIc5iZUalgfz8IAiCKgKm6Bgad1rU71tQAy5ZJ56dOZQdgIiJSPWZmVEqjERDs34a6GYsFmDFDOlksHTQ6IiKizsNgRsXYBZiIiIjBjKqxCzARERGDGVVjZoaIiIjBjKqFBjAzQ0RExGBGxZiZISIiYjCjanIwY2IwQ0REXRj7zKhYmwqA9Xrgyy/rzhMREakcgxkVa9M0k58fMGVKB42IiIio83GaScW4NJuIiIiZGVULbUtmxmoFVq2Szt94o5SpISIiUjF+k6lYm6aZzGbg1lul8+XlDGaIiEj1WjXNlJ2d3VHjoDbgNBMREVErMzMpKSmIiIjA4MGDMWTIEOVksVjwzjvvYOnSpR01TnKCfWaIiIhaGcxkZGRg79692LdvH/bu3Ysvv/wSOTk5AIDQ0NAOGSA1rS4zw2CGiIi6rlZnZlJSUnDDDTcol23duhXTpk3D3//+d3ePjVogZ2YsNjuqa2ww6LQeHhEREVHna/fS7NGjR+Ptt9/G66+/7o7xUCsE+/tBEKTzJtbNEBFRF9WqYMZisTi9vFevXjh06JBbBkSu02gEBPuzboaIiLq2Vk0zBQcHo1+/fhg6dCiGDBmCoUOHIiEhAe+++y4mTpzYUWOkZoQY/FBmtroezPj7A4sX150nIiJSuVYFM7/99hv279+P/fv3Y9myZXjqqadQXV0NALjqqqvw/PPPY+DAgRg4cCD69OnTIQMmRyEGHVBa7frybJ0OmD69Q8dERETUmVoVzIwbNw7jxo1T/rbb7UhPT8e+ffuwb98+7NixAx9//DEKCgpgs9ncPlhqjMuziYioq2tX+1eNRoO+ffuib9++uP3225XL8/Pz2z0wck1dMONiZsZqBX76STo/aRI7ABMRkep1yDdZbGxsRzwsOdHqXjNmM3DttdJ5bmdAREQ+gLtmq5ycmTFxmomIiLooBjMqx/2ZiIioq2Mwo3IsACYioq6OwYzKhba2AJiIiMjHMJhROW42SUREXR2DGZULDZALgJmZISKironrclXOGCBlZkqrXAxm/P2Bf/2r7jwREZHKMZhROSWYqWzFdgazZnXgiIiIiDoXp5lULrQ2mCkzW2G3ix4eDRERUedjMKNycmZGFF0sArbZgPXrpRP3zyIiIh/AaSaV0/tpYdBpUF1jR2lVDYyBuubvUF0NXHaZdL68HAgK6vhBEhERdSBmZnxAq4uAiYiIfAiDGR8QFiCtSiqpsnh4JERERJ2PwYwPYGaGiIi6MgYzPiCUwQwREXVhDGZ8ADMzRETUlTGY8QEMZoiIqCvj0mwfEFa7HNvkSjCj0wELFtSdJyIiUjkGMz5AzsyUuLKlgb8/MHduB4+IiIio83CayQdwmomIiLoyZmZ8QKuCGZsN2LNHOj9sGKDVduDIiIiIOh6DGR/QqqXZ1dXAhRdK57mdARER+QDVTTP94x//gCAImD17tqeH4jU4zURERF2ZqoKZnTt34sMPP8SgQYM8PRSvIq9mKqu2wmYXPTwaIiKizqWaYKa8vBxTp07Fxx9/jPDwcE8Px6vImRnAxeXZREREPkQ1wcysWbNwzTXXYOLEiZ4eitfRaTUI9JcKeTnVREREXY0qCoC/+OIL7NmzBzt37nTp9mazGWazWfnbZDJ11NC8hjFAh0qLjcEMERF1OV6fmcnOzsYjjzyCZcuWwWAwuHSf+fPnw2g0KqekpKQOHqXnsQiYiIi6Kq/PzOzevRsFBQUYNmyYcpnNZsPGjRvxr3/9C2azGdoGvVKeeuopzJkzR/nbZDL5fEDj8vJsnQ6YN6/uPBERkcp5fTAzYcIEHDhwwOGyGTNmoE+fPnjiiScaBTIAoNfrodfrO2uIXiHM1WDG3x944YWOHxAREVEn8fpgJiQkBAMGDHC4LCgoCJGRkY0u78o4zURERF2V1wcz5BqXgxm7HThyRDrfty+g8fqyKSIiomapMphZv369p4fgdZRgpqWds6uqADmjxe0MiIjIB/BnuY8wBnKaiYi8n51dyqkDMJjxEayZISJv98GGkxj84s84nOP7vb+oczGY8REMZojI2/12pABlZit+PZLv6aGQj2Ew4yMYzBCRtyuskDqzH8llZobci8GMj2AwQ0TerrBMCmYOM5ghN2Mw4yPkYKbcbIXVZvfwaIiIHFmsdpiqrQCArKJKlJutHh4R+RJVLs2mxuTtDADAVG1FRJC/8xvqdMBjj9WdJyLqBOcrLA5/H801YURqhIdGQ76GwYyP0Gk1CNb7odxsRWlVTdPBjL8/8M9/du7giKjLKyw3O/x9hMEMuRGnmXwI62aIyFs1DGZYN0PuxGDGh8hTTSWVlqZvZLcDmZnSyc7aGiLqHEXl0ueSIEh/H84t8+BoyNcwmPEhxgBp1rDZzExVFZCWJp2qqjppZETU1RXVLsse1C0MAJCeZ4KN3YDJTRjM+BB5msnEaSYi8jKFtZmZ4cnhCNBpUV1jR0ZhhYdHRb6CwYwPYc0MEXkruWYmJlSPPvEhAFg3Q+7DYMaHhAVKK5gYzBCRt5FrZiKD/NE3PhQAOwGT+zCY8SFGpQCYwQwReRc5MxMVrEe/2mCGG06Su7DPjA8J5TQTEXkpOTMTFaxXPquYmSF3YTDjQ1gzQ0TeSBRFZTVTZLA/jAE6CAJQUGZGYbkZUcF6D4+Q1I7BjA+JqK2ZadicyoGfH/DAA3XniYg6mKnKihqbtAw7IsgfBp0WqZFByCiswIEzpbisT4yHR0hqx5oZH5ISGQgAOH2+sun+DXo98N570knPX0NE1PEKa7MyIXo/GHRaAMDwlHAAwLPfHMSZ4kqPjY18A4MZH5IQFgB/Pw1qbCLOFrMhHhF5B6VeJqTuB9QTV/VB96ggnC2pwh2fbEeBqdpTw6MGamx27M8uwZ7TxdhzuhjZ570/2OQ8gw/RagSkRQYhPb8MpwrLkVybqXEgikBhoXQ+KqqutzgRUQcpqp36jqy3AW50iB7L7h2FKR9sRWZRJaZ+sh1f3T9Gqf0jz3l0xT58/0euw2XfPTgOA7sZPTSiljEz42PSooIAAKfONdFZs7ISiImRTpXeH20TkfrJdXyRwf4Ol8cbA7D8rxchNlSP4wXl+G5/jieG55TVZkeVxdbsqbrG5ulhdogDZ0sBALGhegT5S9OCu7LOe3JILWJmxsekRUvBDNuEE5G3KKy3LLuh5MhAXNEvFp9tO+01U02nzpXjhvd+h6na2uJtp49JxQvX9++EUXWec2VS8PnlzNH4fEc2PthwEple/p3CzIyPkTMzDGaIyFvULct2vuggvHYlZrGXNPzceqrIpUAGAH46lNfBo+lcFWYrKi1SxikqWI+0KKlcIaPIuzP5zMz4mB7MzJAXEkURNrsIPy1/P3VFhWVyZsbf6fVhSjBj6bQxNSe/VMoQ/XlEEuZd38/pbQrLLLjkn+uQb6qG1Wb3mde2nJUJ8tciSO+H1EjpO4WZGepUaVHBAICzJVU+O59L6vPoin0Y/vKvqlgVQe4nZ2aaao4XHuhdW7Hkm6TxJoYHINDfz+mpW3gAdFoBdhHIL2umt5fKnKutb4quXXkmZ/vPFFfCYrV7bFwtYTDjY8IDdcpqAGZnyBuUVFrw3R+5KK2qwYqd2Z4eDnlA/U0mnQn3ssxMXm3tTlyoocnbaDQC4ozS9bklvtMKQ87MyMFMdIhUBGwXgWwv7gfEYMbHCIKA7pxqIi+y4dg5pYnjqr1nYW+qoSP5rHPlzdfMhHldZkYKZmKNTQczgLQaC5Ay4b5CLsKWgxlBEJBSO9WU0dQqWS/AYMYHNVsE7OcHTJsmnbidAXWCtUcKlPNnS6qwM9O7l3iSe5mtNpTVFtNGt1gA7B2ZGSWYCW2+S3pimBTM5JZ6xyosd1Cmmer9W8nfKZlF3hvM8NvMB3VvrteMXg8sWdK5A6Iuq8Zmx/p0KZjpExeCo3llWLX3LEZ1j/TwyKiznK+QAhQ/jYDQAOdfOXIwU2mxwWy1Qe+n7bTxNVRdY1NWVTU3zQQA8bWZmxwfysw0nGYCgFR5RZMXZ/uZmfFBchHwqcJyD4+EurrdWcUwVVsRHqjDs9dIq0L+dyCXxeldiLySKTLYH0ITHcdDDH7Q1F7l6ammgtriX72fpsVuxAm1mZmcEh/KzDgLZiK9PzPDYMYHNVszI4pARYV0Elm7QB3rt6NSVuayC2IwpkckEowGlFVblcvJ9xW2sJIJkIppvWV5dp4yxWRoMviSJYT5YGamvHEwI3+nZBayAJg6kRxFl1TWoLiiwQdDZSUQHCyduJ0BdbBfj+QDACb0jYVGI2Dy0EQAwNd7znpyWNSJlJVMzQQzQF0RcHGFZzMz+S6sZJIlKDUzPhTMyJmZ4Lrjl79Tckq9t+UHa2Z8UIC/FglGA3JKq3GqsALDm1gOSdSRMgorcOpcBfw0Ai7uHQUAuGloIhatP4n16QW47aOtje4TFazH/JsGIsTAzQZ9hbzJZFQLn0NS3UwFSjycmXF1JRNQF8wUV9agymJDgL/nan3cwW4Xla0nYuoVP0cE+SPE4IeyaiuyiipxQVyIp4bYJGZmfBT3aCJPW1ublbkwLQKhtcFJr9gQDE0Og9UuYtup841O3/+Rix8P+FZ7+K5O3mQyKqT5zIzcOK89Wxqcr7DgvXUn2pUpySuVMzPNjxcAQg06BOulnECOD2RniistsNlFCIIUwMgEQfD6rXKYmfFR3aOC8fuJIpw6xyJg6jx//+4wPt9xGiJE1NikmqzL+8Q43Obju0ZgR8Z52BvUbH23Pwc/HcrH4VxTp42XOl5LDfNk7qiZ+XRrJhb+ehw5JVV45caBbXoMuZtvrAvTTIBUN3Msvxw5JVXoER3cpuf0FnK9TESgP3QNtmdIjQzCH2dKvbYImMGMj/L2KJp808pd2aiqN6cerPfDnwbGO9wmKljf6DIAsFjtDGZ8kFxQ21LNTN2WBm0PZo7llwEADua0/TUk78vkajATbwzAsfxy5PrAiiZnK5lkqVHevUcTgxkfxWkm6mzlZivKzFJztF/nXAKDTovwQH8E6V37mOmXEAoAOJJjgiiKLa4kIe9ntdmxP7sEgNRnqDlhbtg5W+6tdTy/DHa7CI2m9a8hZSsDF2pmgLq6GV/oAtxcMJPm5b1mGMz4qOQI6YV3prjKo18MeaXV2Ff7YdYW3cIDMCDR6L4BdXG/HM5HkF6LMT2i3P7Ycq1BiN4PPWNaXyDYIzoY/loNysxWnCmuQlLta5jU61COCRUWG0INfugbH9rsbeXGeW3NzNjtovJFW2mxIbu4UmnD7ypRFFu1mgkAEuT9mXygZqZuJZOTzIyX95phMOOj5Dbb5WYrSqtqlF890GqBW26pO9/Bpny4Bdnn2/cm/3XOJW36ciRHheVm3PfZbggA1j12qduDhdasAnFGp9Wgd1wwDp414VCOicFMrV8P5yM6RI/BSWGeHkqrbc8oAiAVgWtbyJK0twD4bEkVzPV2dT6SW9bqYKa0qkZ5DGfZCWd8qXFe85kZ6f9lvsmMCrPV5YxrZ/Gu0ZDbGHRaRIfoca7MjDPFVXXBjMEArFzZKWOwWO1KIDM0OQzaVmaH0vPLUFZtxbH8cgYzbnCmuErZ8PH99Scx/6a2FUg2pW4VSNuCGQDoFx+Kg2dNOJxrwlUD4tw1NNXKLKzAvf+3C8YAHXY9MxF+WnUtQN12StqH6yIXtq9obwHwyQaLHdLzylr9GpKnmMIDdTDoXPuxFy83zvOFzIyThnmysEB/hAfqUFxZg0+3ZiHOqEdkkB7jeka1aTrP3RjM+LBu4QG1wUylR6ZqSqukX1iCAPz3vjGtfsHf/9lu/HgwT9nFldonr95meF/tzsaDl/dUMnhueXxT6wonnZGnIg63o4CzKRmFFXj712NKgbJOq8G1gxIwqX+s19bnHMwphShKDTAPnC3F0ORwTw/JZTa7iJ0ZUjAzKq3lYCY8qH07Zzfci+5oXutfQ3mtLP4F6rLgOSWendJ3h+YyMwDQPToYu7OK8dqao8pli6ePxGUNVix6AoMZH9YtPBB7T5fgTLFnfjGUVkm/sEINujZF7jG1b6iC2jcYtU9BWV0wU2MT8f66E21evuqMUmtgdC0970y/2mDmSAesaHpv3Ql8sy/H4bLv/8hFv/hQPDyhF3rGSMtqDToNuoW3b4rLbLXhrV+OK4G4IAiYMqKbSxmK+o7n12UbtpwsUlUwczjHhDKzFSF6P6W4uzn1a2baUrwr70U3qJsRf5wpRXpeWavHLO/L1JpgRi4Urq6xo6SyBuEqblLaXM0MADw6sTc+2XwKNruIY/llyDeZcaKgnMEMdaxu4dIvBodgpqJC2soAAMrLgaDWzSm3hvwLS25T3loxtR8o5xjMuIX8q7N/QigO5Zjw5a5sPHCZ+7IzuW6YZupb+6V3tqQKJZWWuulRJ6osNmzLKIK1tp9NkF6Li9IinX4JiqKI308UAgBmju+O5IhAnD5fic+2ZuFwrgn3fbbb4fbPXtMXf724e5uP47cjBfhgw0mHy/ZlF2Pt/7u0VY9zvKDuC3nLyULMuqxnm8fU2eR6mZEu1MsAdZ8TdhEoq7bC2MrPDTkz86eB8fjjTCkyiipa3ZU3r5XFvwCg99MiKliPwnIzzpZUqTqYKWghMzOuVxTG9ZIWD7z0/WH8e3OG0hTR0xjM+LC6YMYzezApwUwLO882Rf51wMyMe+TX/ur808B4hBp02HqqCIvWn8DLN7gnO5PvhmmmUIMOSREByD5fhcO5pmZXXT3+3z/w3X7HTMtLk/vjztGpjW57qrACuaXV8PfT4NGJvZV6iJmX9MAnm07h6z1nUW21wVxjR1WNDXtOF7f5GOTnA4AhSWGY1D8Or605ipPnKnC+wuLQWbUlx+plZnZlFqO6xuZyLYenbTslBTMXdY9w6fZ6Py0C/bWotNhQXGlpczBzYVoEIoL8cb7CguMFZRjULczlx8hrYxF7QpgBheVm5JZWq3b1pdlqU0oDXCl+ljcOPeclwYy6qsmoVeRUueemmaQ3Rmhbg5lQBjPuVH/J6SMTewEAvtx1BjU2e3N3c5mc+Yk3ti/T08+Fupk/zpTgu/05EARgcFIYetT2Vfq/bVkQnewGL2dlhieHOwQDEUH+ePyqPtj29ATse/5K/HPKIABQ9qdpK7mx2OV9YnD/pT2U8e1tRZBksdqVxwnQaWG22rH3dEm7xtVZbHYRO1pRLyMLb2MRcIXZqgQiPaKClZ42R1s51VTQhswMACQY6+pm1Ep+zftrNTC68JkdFezvcD9PYzDjw+pPMzn7gO9oJVXyNFPb0q5yzQynmdyjfoHuhakR8PfTwGK1OxQGt5XVZlfSzbHtqJkBgH7x0i/b5joB//OndADADUMS8e2ssfj6gbEw6DQ4ll+OvU76GsnBjJwib0pkkDT29qbOs85L2dCUSOkHxfAUqdZld5brwUxmUQWsdhHBej9c0S8WALD1ZGG7xiV7dMU+3PbR1g7bAflIrgmmaiuC9X7o70K9jCwssG1FwHJ/mcggfxgDdcpGiEdzWxfM1L1HWvca9oUVTfWLf10pYpb32ir0ks9nBjM+rGGvmc5WWvvrqs3TTLVvlqIKM6xuyh50ZfULdDUaQWn25Y7M3blyM+wi4KcREBXUzmAmofnMzJaThdh0vBA6rYBHJ/YGABgDdPjTAGmLhC93Zjvc3mYXseWkNOUxtmfzwUx0iBR4F7Xz12ZWbWMxuc9JW4IZuTV/z5hgjKsd9++1x9Eep86VY9Xes9h26jy+3nO23Y/nzPbarMyI1PBWLSdva2ZGXpbdvTYDJmdm0vNbV0ieV9r6AmCg/oom9a68lIOZljYElUV72TQTa2Z8mEFXV5jm0Gumk9RlZtoWzEQG6aERpILA8xUWpSCYWq/SYkVZtbTVgPxBnRgegMyiSrekxuXi35gQfbt7TvSNl76IThSUw2y1Qe9XNy0kiiIWrJGyMrdfmIzkyLpVR7eOTMLXe8/iu/05eO7afkpTrwNnS1FWbUWIwQ8DW6hnkDMzpVU1sFjt8Pdr/e+9SotVqU9KbZCZ2X+mBDU2e6NN/JyRVzL1jg3G6B7SVM3+7BKUm63KTs1tsfZIgXL+P79n4LaRSW7vE7JdqZdp3eqtsDY2zpPrZbpHSYsb+sRJAXFrMjM1NjuKKqR/N1e3MpDJU6vrjxZg8nu/N7r+hiEJmDE2rVWP2dlaWsnUkPxj83xF21afuRuDGR/XLTygNpjp/F4zcqrYlflXZ7QaAVHBehSUmVFQZmYw0w7yl2ugv1b5Ikx0454yyuZ8bez+W19iWABCDX4wVVvxxFd/IMRQ9/oprarBvuwSBOi0ePByx5U9o9IikBoZiMyiSvzvj1zcOjIJQN0U0+jukS2uqjEG6KDVCLDZRRRXWpz+Qj9RUI7pi3cor2+tRsBjV/ZWCo9P104xGQN0yg+I7lHByjEdzS3DwG4tvxfllUy9YkKQFBGoFEbvzDyPyy5o+1LYX47kOxzLhuPn2vV4zuw/UwIAGJHSuqXkbd3SQC64ljMzvWNDIAhAUYUF58rMLhW0niszQxQBnVZARCt/+PVLCIUgAGVmq7IXVX1Hck24/cJkry7ebqnHTENyIbv8XmlpI9GOxmkmH9doebZWC/zpT9Kpg7czkKe22hrMAHVvrPo9Uqj16nfnlefDE8OkrMFZN0wztWVJa1MEQVBa93+zLwf/ty1LOa2uXb1097hUxIQYGt1PDmBW7KqbanK1XgYANBpB+ZBuqm7mu/05OFNchXKzVZnCXbb9tHJ9VpEUzKTWyxppNAKGKVNN51scB1C3kqlXrJRtGFu7smtrO6aaiissylTXpP5SHc5/Nme0+fGcKSirRr7JDI0Al/rL1Fe3pUErg5naaaYe0dL/qwB/rbKXUEvN82psdlSYrcq/W0yIodVZhrSoIPz4yMX497QRjU7RIXpYVFC8fa5ceg+7GszotBrl38sbioCZmfFxjVY0GQzA//7Xqsew2UWUVTtP+zbXEK+9BcCANG1xCF2jCDijsALTF+/AlOHd8ODlvdz62M6WTSeGN5+ZySiswHvrTjTZfKxPXAj+cfMgaDVCq3cabsnLNwzA6n05sNobF66HGPxwx0UpTu93y7BueOPnY9idVYzdWefRL96IXZnSl3dL9TKyyCB/nCszN/kBLWcd5lzRGxf3isKN72/BsfwypadJw3oZ2fDkcKxPP4fdp0swfWzzY6i/kql3rDTtNrpHJL7YmY1l27KwPr2g0X1GpEbg5ckDmv0iXn+sADa7iAtiQ/DsNf3wy+F8bDpeiKN5JmVqpr0Oni0FIAUWgf6t+4ppy87ZdrtYN80UXff//ILYEGQUVmDx75nYmVkMvZ8Gtwzv5vAeOHi2FH/+cCsqLHWF0K0t/pX1iQt1+v9wTI8cfLsvB1tPFirThd5EFEXYxbrsravBDCAtzy6urEFhuRkXwLNbzjCY8XHt7TVTXWPDVQs3IrPI+f0HJIZi9axxTj9AlQLgNtbMAFB+fcudOX3Z++tOIKuoEt/tz+3AYKbug6qpaaackios/PUY/rvnrLKXkzMHzpbi1pFJGJkaoUwzuSMzA0iBwEMTWv//ICbUgMsuiMGvR/Jx86Kt6BMXAovNjnijAd2jXGsQGR2ix9G8MhQ5ycyIoqhMI4zvHY1B3YzKHmiHc00YnhKuvFdSIh27CMt1M3tcKAKuv5IpvjZAHNczCoH+WlRYbA79Z2TH8svRNy7EaZ8d2a+19TIT+8UgKSIQVw2Iww8H8rBo/UnMuUIqpo4K1rdrE8EDZ6RMSEv1Sc7UbWng+i/9PFM1qmps8NMIDpuT9k8IxZpDefjtaAF+Oyodd/b5Svzj5kHKbTYdL3QIZAQBuLK/e/cEG9MjEt/uy8GWk0WY49ZHbpv5PxzB8u2nUWO3w2oTG/1gcLVmBpBeK8cLyr2icR6DGR/ntAtwK/x+orDJQAYADp414Vy52WltgZKZccs0k+ffLB3pXJkZ39a22s/vgCk1Z83AutXLzMgFfFabHbd/vE1JuU/oE4NbRybBv0HB6sebTmHLySJsP1WEkakRbs/MtMfLNwwAIOK3owVKn5ExPaJc3jMnMqjpFU3Z56tQXFkDf60GfeJDIAgCBiUasfZoAQ6cKcHwlPAmMzODk8KgEaT/37mlVc3246m/kkked2SwHr/MGa88fn3bThbhnd9O4NUfjmJ87xiHwmiZxWrHxvRzAICJfaUppnvGpeGHA3n4dl+O8vozBuiw/rFL29zJ9kBtZqYtNXpKZqbC9cyMnJVJjgx0KKy+c3QKqmpsKKu24kxxJdaln2u0GaVc/D5zfHc8OrE3BAEOBefuIDd+3Jdd4vHdpistVvzn9wzU2Jz/SIkM8sew5DCXHy/Ki9pnMJjxcfWnmURRhFBZCcTUFvsVFLS4nYH8i+aOi5Lx4vUDHK4b+4/fkGeqRl5pdaNgxm4XYXJDzUxMqPe8WTrSZ9uyYKldfl5SWeP2Tq/KnjP16kzijAYIgvQlV1hhRkyIARmFFcgqqkSATovP/jpKySY0dPp8pRTMZJzHg2jbBn0dJc5owCfTRuJMcSU+33Eae0+X4G+XuL41gVzIWFjR+DW3r3aKqW9CqPKlN6A2mPmj9ks8s7BxzQwABOn90CcuFIdzTdiTVYJrBjUdzNRfyVRfYliA0+0nLkqLxPaM89iecR6PfbUfX9x7UaNs6Y6M8ygzWxEVrMfg2q64w5LDccOQBPxyWCoKrqqRusBuO1WEqwfGNzm+5sjTTK4UOTfkagHwoZxSLFiTjuoam5IVkFcyycICpYaIgNSscF36uUb1YXJWMiUiqMOKc5MiAtEtPABniquwK6sY43tHd8jzuGJ7xnnU2EQkhgVgxcyL4KfRQKsR4KcRoNUKCNRpW7WU3psa5zGY8XHyr2+5UDEMACpdm3ISRVEJZib0jW20EiTWaJCCGVM1Bje4b5nZCjl72dYOwED9zSZbl604kmvCf3efQffoYAxNDkPv2BCX9ofxhOoaGz7bluVw2bkys0PKvL2cZU50Wg1iQ6R/w7PFVYgJMSiZjD7xIU0GMgAwqrZF/e6sYtTY7G4tAHaXbuGBmDupT6vvFyl/QJc1/oCWp5iG1PuiHlR7/uDZUpitNuTWNk5rmJkBpKmmw7km/HQoT3me+vrGh8IYoHNYyeQKjUbAP28ZjKve3ogdGeex8NdjuH5IIhLDApS9iX6tXcU0oU+MEugIgoCFtw1VHueZVQewbPtp7Dld3KZg5lyZGXmmaghCXSfn1gh3cWn2K/87ovQOkg1Jajp4kuvD8kzVDkvj5cxMQljHvm5Hd4/Eyt1nsOVkoUeDmc3HpWL4i3tFtXszVaDelgZe8GPT64OZ+fPn4+uvv8bRo0cREBCAMWPG4LXXXsMFF1zg6aGpQqNeM2Gu/5MfzStDbmk1DDoNRjvpFxEXqsd+1NVj1Fda+2EUoNO26xdPW6eZXvnfEWw+UdctNTxQh+X3XoS+bfiA7Wir9+WgqMKCBKO00uhsSRXyTdXuDWaayJwkhgdIwUxJFYYmhysrP1oqBu0dE4KwQB1KKmuw5WQRqmukrJI3TDO1l9z0r8hJZkYOZuTVVkBdbciJgnIcyyuHXZSWwEc5CVaGp4Qrq7JWN9hXCgBCDX548uq+StF1rwaZmeYkRwbiqT/1xXPfHMQ7v53AO7+dAAD4+2kgAErmb0LfppdhD0sOrw1mSlx+3vrkrEz3qKA2TafI00xVNbYms5NHck3YcrIIWo2ABTcPgkEn7ek0pmfTxbVRQXqHjtfye0vO1Mg/+jrKmJ5SMNOelWjusOm4NM3oyso+V8j1NayZccGGDRswa9YsjBw5ElarFU8//TSuvPJKHD58GEEduOOzL3HoNRPm+pe5nJUZ2yPK6YeK/CvcWTv8kqr2F/8CdQXAUg8I0aW6B7u9rkhzaHIYjuaWobiyBuvTz3VKMJNvqm5y9Zcz/65dGjttTCp+PZJfG8y478NBFEUls9VwpUZiWAB2ZxUrH+pykzG5cV1TNBoBI1Mj8MvhfHy7T+oiawzQeXUfDVfJGZOGNTNWmx0Hc6Qv6/rBTEyoAbGheuSbzPjhYC4AKSvj7LU6sV8sLu4V5bRRYblZarb39KoDymXySiZXTb0wGfml1fj5cB7OFlehwmKDxVrXPTvBaMDFvZrODMjLxw/UZplaWz8i18u0pfgXkII5uc9PaVWN09fTkt8zAUhLy28e3s2lx9VoBCSGBSCjsAJniquQFBEIU3UNysxSI8kEN+0c35TR3aXg4eDZUpRW1bRr6r2t8k3VOJZfDkGoW+bfXvKPTQYzLlizZo3D30uWLEFMTAx2796NSy65xEOjUpdu4QHYl10iFQGntT6YubyJX3JyMWmes8yMG+plgLo3i9lqh6na6tLjZRRVoMxshUGnwcqZo/HWr8fw3rqTnbIJ3LqjBZixZGer7xfor8VtFyYrdRfOsl1tdb7CohT8NezN0nB5tjLN5MIy3VFpUjDz8yFp+iLeB7IyQF3qvOFqpmP55aiusSPE4Ie0BlNIAxPDkG/Kx/d/SNmWhvUysmC9H/7vnlFOr7Pa7Fi6NQtv/pyOCosNIfVWMrlKoxHw2KQL8NikCyCKIkxVVpSZ6wLrqGB9swFnamSgsuP0oRwThiU33/Rud1Yx9mQVY/rYVOi0mnYV/wLStFdYgA5FFRanTQvPV1jwTW3w3NqOunIwI7/W5c+D8EBdq5eQt1Zc7Wq6U4UV2JFxXtlrqzPJU0wDE41tLu5uKMqLMjOqa5pXWiq9WSIiXNtWntq2e/b5Cgv21O7w21R3UDkz4+yLt73df2UGnRYhBumDxtV52T9qizT7Jxjhp9Uov7oaBjNrDuZhwhvrldS4O6w9Kn2xG3QahAXqXDpFBvljzhW9YQzQKQW67lzRJGd5ooL9G7XnV5ZnF1ehtKpG+aCXN+prjtyqvtzsuE2C2ik1MxUWhw1a5f4yg7uFNSquletmss83XS/TEj+tBveMS8Ov/288po5KxnPX9XN5BZYzgiDAGKhDt/BA5dRS5kwQBAytzTq1tITcZhfxwLLdeOWHI/hX7ZTWwXZmZoB6Wxo4WdH0+Y7TMFvtGJAY2uruwvVf6/X/29FZGZncY2aLmzYLbS152n2ci/2WXBFVby8zezNtHDqD12dm6rPb7Zg9ezbGjh2LAQMGNHk7s9kMs7nui89kat1mY76mLb1mNhwrgChKBYlNvdmbn2Zq375M9cWE6FFWbUVBWTV6xrRcQ/DHGccP1Kb6qfx3zxmcPFeBX4/ku22rhwNnpdfaglsG4/rBCa2+v7x6y519deRgs2FWBnDMzMh1GolhAS4FoX3jQxGi91NS9d5U/Nse8v5MFqsdZWYrQmu3U6irl2n8Wmm4cqdhj5nWiDcG4JUbB7b5/u01LCUca48WtNixduPxc0qg/P76ExjVPQK5pVLxb/92vJ+kFU0VjVY01djs+HRrJgDg7rFprQ706loRSJ+D8o8bZ6vDOsKYHlFYtv00Fv+eicW1U2X1XdwrCktnXNghexyJoohNx13vhO0q+b1irZ0WdFfGpy1UFczMmjULBw8exObNm5u93fz58/Hiiy920qi8n/wm/v1EEd5aexyzxl0s9Q3RNJ2Ykzeju7xP0/Pr8jSTs/qOuh2z2//ijgkx4OS5ilZkZuS6BsdgpmFmJrt2D53zFe5ZVlhjs+NIrhTMDGrjh7lc0+LOaabmesB0qxfo1RX/ulanodUIGJEajnW1vUvcsS+TNwjw1yKotjldUblFCWb2ycFM7bLm+hpmItoTzHiaPLUkZ2ab8tWuMwAAf60GFpsd93+2B4DU2r89G2HKRcCv/ngEH2w8pVxebbEh32RGVLAe1wxq/UqrhlOqZ0o6NzMzrleUUlvlzKbjhfjfgVxc14YfQS05mleGwnIzAnTaZlcptpa/nwbGAB1Kq6QuwAxmXPDggw/i+++/x8aNG9GtW/NFX0899RTmzKnrtWgymZCUlNTRQ/Raw1LCcUFsCNLzy/D2lrNYNP4p9I4NhvDJ7ibvk17btOvyPk3P7cq/xOU9aup/gJW6MTMT3YrGTFabHYdqizQH1X7pxNd+WJmqrSirrkGIQQdRFJUNAYvcFMwczy+HxSrVVLT1y0yZZnIxmDl1rhw5Jc5v2yc+BFHBeqdbGcjkD/iyait2ZJxX7ueqUd0jlWDGVzIzgNRrpuJ8JYrKzUiLCkKlxao0shtSr/hXFhWsR4LRgJzaLGVqG6aZvMXgJCO0GgG5pdVNNvcrqbQovWkW3TEMj3yxT3nPt2eKCZBWcP16JB/Z56uUabv6ZoxNbVNjO/lHjTzdLr9vOnolk8wYoMPmJy5X+m/Vt2RLJt797QQW/noMfxoY7/Y2EnK9zKjuEW5vChgV7I/SqhqcKzejVysL1t3J64MZURTx0EMPYdWqVVi/fj3S0lou+tLr9dDrPbuDpzcJNejwwyMX45fD+fhw40nsPV2Cg2dbnnpLDAtw+sEtC9L7KdMMeaWOU0ByzUx7eszIYlqxPPt4QW2Rpr6uSDNY76f8esgtrUaIQSowrKxtY17spmBGrhcYkGBsc62DvDO4K8d6oqAcVy3c6HT/IkBaufLTo5c43cpAFujvh/BAHYora7DhmBSUtGaPnlFpdbVrcUbfec9FBfvj9PlKpRnYwbMm2EWpyLmp3dsHdjMip7Qa/n4aVQd2gf5+6BMXgkM5TTf3W70/BxabHX3jQzGhbyyevLoPnv3mIID2BzOPTuyNMT0iYa6xN7ou0F+LUU7aRLhCDtxzS6pht4s4Wzvt3lmZGUDq7eRsd+m/XdId/7ctCyfPVeDbfWdx07BuKK2qwdNfH1BqtZrjpxGUJeo6rQYNP35O1nZJdme9jCwqWI+T5yo83jjP64OZWbNmYfny5fj2228REhKCvLw8AIDRaERAQOe9CNVOqxFw1YA4TOofi0M50hYELemfENriL4RYowFlBeXINzUIZtxZM6PUkbScrZCLfwckGh3mnhPCApQC196xIUpWBnDfNNMfZ6XnbkvnU5kccJRVW1FpsTa7ymLV3jOw2kVEBPkrAZ8sp6QKOaXVeP2ndCWt3dQXbEJYAIora1BWLdW+tLQsu74BiUYE6/1QbrYiyQ1NuLyF/IUj95rZ38wUk2xQtzD8dCgfyRGBHVL30JmGJYfjUI4Ju7OKnU7prKydYppSuzT6Lxcm4+fD+dhyohCXtLMpnL+fptnl420VF2qAViPAYrPjXLlZycx0Vs1Mc0IMOvztku5YsCYdb689jou6R+LuJTuV1YXuoBGAy/s03WOorbxlSwOvD2YWLVoEALj00ksdLl+8eDGmT5/e+QNSOaGyEgOG9JT+yMxscTuDlsSFGnCioLxREbDcNM8dNTPKNJMLAZhcLzOoQUCRGGbAkVyTUjeTXS+Ycdc0k1z8255i4mC9HwL9tai02FBgMiM1yvlbVBRFpenai9f3bzTPvvl4Ie7493Z8ui1LKeZtarVRYlgADuVIY/f307RqikSn1eDdvwxF9vlKj6aY3S2qQRdgeRuDwc1kKi+7IAYLfz2GSzrgi7izDUsJw/9ty3JaN3M0z4QDZ0uh0wq4YWgiAGlJ+H+mjYCp2ooID9ZNNMdPK2XMzpZUIaOwQlkx2JmZmeZMG52Kf2/KQFZRJa58ayPKzVZEh+ix4JZByjYPTbHa7Ki02FBpsaHG1jijBUh1XN2jXW/C6CpvaZzn9cFM/aWR5CaF7lsaKH9BNuw1486amdbsnF0XzIQ5XN5wefbpeptnFtcuwW3PMtj6xb/tSbMLgoDYUGmPpHxTNVKb2Ol5b3YJss9XIdBfq2waWN+4XlG4cWgiVu09q0z5NRnM1KsZ6B0b3Kq9WYCml+6rWWRQE5mZZlrm90sIxf55VyLABxoHDk+Wpg8P5ZTimXpN/AAor/MJfWIdAhc/rcZrAxlZYlgAzpZUYXdWMURRCt6ddWr2hCC9H+6/tAde/t8RlJutSI4IxP/dc2Gblvl3JqVxHjMzpGZynUTDglW5A7A7Ol26WjNjttqUFTkNMzN1wYw0zvrTTFa71FzM2I7ASyn+1fshpZ3bEMSE6Gt/OTZ9vKtrdzi+sl+ssvdOQ89e0xfr0guUYKaprQbqp9lbUy/jy+p3AZa3AhGElgPVjm6+1lmSIgIQbzQgt7Qay7afdnqbW0e61n3Xm3QLD8COTCjF7olhAe36EeNud1yUgp8P58NPI2Dhn4c0WZ/lTeo2m2QwQyoW20SvGXc1zQPqIv/SqppmW6wfzS1DjU1EeKCu0QqFhAa9ZuoHM4D0C7w9wczBep1P21svoRQBN1EjZLOL+P4PqW3+9UOaXsYZGazH01f3xeP//QMBOq2yiV9D9f9fubos29fV72wq12H1jA5GiKHz29B7giAI+OjOEViXLvWbaighzKDKjJychdxd2xDQG+pl6jPotPhy5mhPD6NV6t4rLAAmFYt10gW4usYGc+1+MO6YZjIG6JRN4s6VmZvc7VX+0hnYLazRr63E2l1xndXMAFIRcPd2lDooe9K0o/hXFhvSfK+ZrSeLUFhuRligDuN6Nj/oKSOkVREJzfwCTQyr+//pjRtxeoKSmamwYF924/2YuoKB3YxueT17Ezl4KVf2ZPL+zIe385YtDRjMULvEOamZketltBqhXc2zZIIgIDpYj7MlVbjhvd/h10SzP3lzx8FOPoDlzExeaTWqa2zIrR1vUkQAss9XtXtFU3v3pKmvLkB0/uGwer+0N83VA+IbbU/QkCAIuPeS7s3epn5mxpVtDLqC+h/QdSuZfOuLvStKbJCxrR/IU9vIq5mKyttfe9geDGaoXeQ6jHNlZlhtdvhpNQ5TTO56YQ9JDsPZkiqXUpmXXtA4WxETIi3LtNpF7DktFf8F+WvRKyak3cFMjc2Ow24o/lXG2kwXYLPVhh8PSu0JJjczxdQa4UH+ePyqC6AVBOVLvKuLrC1kLamswd7aFT1dLTPjixpOKzEz037ye8Vis7e79rA9GMx0NRoNMGJE3fl2igrWQ6sRYLOLKCy3IM5oUPZUCXPjNvdv/3kIZl3aE/YWVrdFBPk7XWqp1QjKssxtp6Tiv6SIQGX1havLs0VRbFRDcCy/zG3Fv0BdZsZZwfPOjGKUVVsRE6LHhanu22z1gUt7uu2xfEFYoD80AmAXpc7R/loNi6N9QMPPhoaZGmo9eTPgsmorzpVXM5ihThIQAOzc6baH02qkKaA8UzXyTNVSMFM7zeTOF7WfVoN+Ce37MpGXZW47VQQASI4IVH5VuJKZ2XDsHO7/bLfSObih/omhbmmW5qwOSSYXLl7UPVL1jdm8mVYjICJIr9QB9EsIbXFKj7yfQadFVHDdv6u3FQCrVXSwtBnwuTILenqoLpzvTmo3eYNBeUWT0mPGjZkZd5BTyvtqdwNOjghUNkZzJZj5+VBek4GMIADXDHLPtI+8FL3SYlMKFWVyE7NhyWFueS5qWv3+I81t60HqIteICULT7QqodeS6GU8WATMzQ+0WF6rHfgAFtR01S924LNud5JSypbZDZnJkIAy1Dc5cmWbKrt2g7oXr+uH6IYkO1+m0gtuW7dbf8yrfVI3g2q6ddruo1G8Mc+POt+RcZL1gprlmeaQuieEB2Jddguhgvds3XeyqvKELMIOZrqayEujXTzp/+DAQ2P4aj7gGvWbkhnlhLbTg7mwN58uTIgJhr92k8XxFy29CeTn3BXGhHd7pNDpUj7JzUjDTozaYOVVYDlO1FQadhkuoO4HcBRhofk8mUpdutZ8DrJdxn/sv7YG7RqegR4z7t0twFYOZrkYUgaysuvNuoEwz1dZ4uLNhnjs1DGaSIwJhqp0SK66oafa+0i67UmYmKaLjPwRjQww4da7CYQuHPVklAIBBiWHQtXLLAWo9OTMTavBr1X5V5N3kDXF7dsA+RV2VO1pStBeDGWq3uAYFq+7cMdud6hf7CYL0t652RVdRC5mZ/LJqWGx2+GkExBs7IZhxsjx7b7Y0xTQ0JazDn5/qes0MTgpjsbUPmTwkEf5+GozpEeXpoZAbMZihdms4zWTy0mAmvl6xX1yoAQadFhG1v76ra+yotFib3FtH3pgyMTwA2k74YnPWOE/OzAxLZr1MZ7h6QBw2HDuHu8eleXoo5Eb+fhpMblDzRurHYIbaTZ5mkr94vXWaKcSgQ6jBD6ZqK5Jq+8EE+WuVrRKKyi0IjHD+lpCLf5Oa2ErB3eT9mfJri6pN1TU4VlAGgMFMZ+keHay6fXKIuipOvFO7yZmZcrMV9366C7ml0he/McC7CoCBurqZ5NpgRhAERAS2vDxbLv7tjHoZoG6a6WRBOex2EfuzSyCK0vPLG28SEZGEwQy1W5DeD/df2gMaAfjlcL6y5YC3TTMBdXUz9TMsES70mqkLZjonMzMsORx6Pw2O5pXh402nOMVERNQMBjNdjSBIS7P79ZPOu8kTV/XBT7MvwcS+UvtHf61Gaf7mTW4dmYQBiaH408A45bL6OyQ3Jbu4NpjppGmmhLAAzLuuPwDgnz+l45t90uaSDGaIiBpjzUxXExgIHDrUIQ/dKzYEn0wbiX3ZJbCLotuayLnTpP5xmNQ/zuEyOTNT3Ewwc7qTMzMAcPuFSfj9ZCH+90cuMgorAABD2fmXiKgRBjPkdmpr/d7SZpPVNTaluDm5E4MZQRAw/6aBOHCmFKfPV7JZHhFREzjNRF1e3WaTznvNnC2RCpqD/LUI7+Q6oFCDDu/ePhShBj/8aWA8m+URETnBzExXU1kJjBwpnd+50y3bGahdRG3b+qYKgOtPMQlurDNy1eCkMOx69gru2kxE1AQGM12NKEp7MsnnqcVppjMeqJdpiIEMEVHT+AlJXV5LS7M7u2EeERG1DoMZ6vKUYKa8iWmmos5tmEdERK3DYIa6PLkAuMxshcVqb3S93GOmM1cyERGR6xjMUJdnDNApm0cWVzbOzniixwwREbmOwQx1eRqNoCy5Lmow1VRaWYOyaisAoFs4p5mIiLwRVzN1NYIApKTUnScAUt1MYbmlURGwPMUUFaxHoD/fLkRE3oifzl1NYCCQmenpUXiduuXZjo3zTnfybtlERNR6DGaIUBfMvP5zOpZsyVQuP1cmBTdclk1E5L0YzBAB6B0bgh8O5CH7fBWyz1c1un5QN6MHRkVERK5gMNPVVFUBl1wind+4EQjg9AkAPHBpT4xIiUBVja3RdUH+WlyYFuGBURERkSsYzHQ1djuwa1fdeQIgbRcwrleUp4dBRERtwKXZREREpGoMZoiIiEjVGMwQERGRqjGYISIiIlVjMENERESqxtVMXVEUV+0QEZHvYDDT1QQFAefOeXoUREREbsNpJiIiIlI1BjNERESkagxmupqqKuDSS6VTVeM9iIiIiNSGNTNdjd0ObNhQd56IiEjlmJkhIiIiVWMwQ0RERKrGYIaIiIhUjcEMERERqRqDGSIiIlI1rmbqigIDPT0CIiIit2Ew09UEBQEVFZ4eBRERkdtwmomIiIhUjcEMERERqRqDma6muhq45hrpVF3t6dEQERG1G2tmuhqbDfjhh7rzREREKsfMDBEREakagxkiIiJSNdUEM++99x5SU1NhMBgwatQo7Nixw9NDIiIiIi+gimBmxYoVmDNnDubNm4c9e/Zg8ODBmDRpEgoKCjw9NCIiIvIwVQQzb775Ju69917MmDED/fr1wwcffIDAwED85z//8fTQiIiIyMO8PpixWCzYvXs3Jk6cqFym0WgwceJEbN261el9zGYzTCaTw4mIiIh8k9cHM4WFhbDZbIiNjXW4PDY2Fnl5eU7vM3/+fBiNRuWUlJTUGUNVh6AgQBSlU1CQp0dDRETUbl4fzLTFU089hdLSUuWUnZ3t6SERERFRB/H6pnlRUVHQarXIz893uDw/Px9xcXFO76PX66HX6ztjeERERORhXp+Z8ff3x/Dhw7F27VrlMrvdjrVr12L06NEeHBkRERF5A6/PzADAnDlzMG3aNIwYMQIXXnghFi5ciIqKCsyYMcPTQyMiIiIPU0Uw8+c//xnnzp3D888/j7y8PAwZMgRr1qxpVBRMREREXY8giqLo6UF0NJPJBKPRiNLSUoSGhnp6OEREROQCV7+/vb5mhoiIiKg5DGaIiIhI1RjMEBERkaoxmCEiIiJVYzBDREREqsZghoiIiFSNwQwRERGpGoMZIiIiUjUGM0RERKRqqtjOoL3kJscmk8nDIyEiIiJXyd/bLW1W0CWCmbKyMgBAUlKSh0dCRERErVVWVgaj0djk9V1ibya73Y6cnByEhIRAEAS3Pa7JZEJSUhKys7O7zJ5PXe2Yu9rxAjzmrnDMXe14ga53zL5yvKIooqysDAkJCdBomq6M6RKZGY1Gg27dunXY44eGhqr6xdIWXe2Yu9rxAjzmrqCrHS/Q9Y7ZF463uYyMjAXAREREpGoMZoiIiEjVGMy0g16vx7x586DX6z09lE7T1Y65qx0vwGPuCrra8QJd75i72vF2iQJgIiIi8l3MzBAREZGqMZghIiIiVWMwQ0RERKrGYIaIiIhUjcFMO7z33ntITU2FwWDAqFGjsGPHDk8PyS3mz5+PkSNHIiQkBDExMbjhhhuQnp7ucJvq6mrMmjULkZGRCA4Oxs0334z8/HwPjdi9/vGPf0AQBMyePVu5zBeP9+zZs7jjjjsQGRmJgIAADBw4ELt27VKuF0URzz//POLj4xEQEICJEyfi+PHjHhxx+9hsNjz33HNIS0tDQEAAevTogZdeeslhzxc1H/PGjRtx3XXXISEhAYIg4JtvvnG43pVjO3/+PKZOnYrQ0FCEhYXhnnvuQXl5eSceRes0d8w1NTV44oknMHDgQAQFBSEhIQF33XUXcnJyHB5DTcfc0r9xfffddx8EQcDChQsdLlfT8bYGg5k2WrFiBebMmYN58+Zhz549GDx4MCZNmoSCggJPD63dNmzYgFmzZmHbtm345ZdfUFNTgyuvvBIVFRXKbR599FF89913WLlyJTZs2ICcnBzcdNNNHhy1e+zcuRMffvghBg0a5HC5rx1vcXExxo4dC51Ohx9//BGHDx/GG2+8gfDwcOU2CxYswDvvvIMPPvgA27dvR1BQECZNmoTq6moPjrztXnvtNSxatAj/+te/cOTIEbz22mtYsGAB3n33XeU2aj7miooKDB48GO+9957T6105tqlTp+LQoUP45Zdf8P3332Pjxo3429/+1lmH0GrNHXNlZSX27NmD5557Dnv27MHXX3+N9PR0XH/99Q63U9Mxt/RvLFu1ahW2bduGhISERtep6XhbRaQ2ufDCC8VZs2Ypf9tsNjEhIUGcP3++B0fVMQoKCkQA4oYNG0RRFMWSkhJRp9OJK1euVG5z5MgREYC4detWTw2z3crKysRevXqJv/zyizh+/HjxkUceEUXRN4/3iSeeEMeNG9fk9Xa7XYyLixP/+c9/KpeVlJSIer1e/PzzzztjiG53zTXXiHfffbfDZTfddJM4depUURR965gBiKtWrVL+duXYDh8+LAIQd+7cqdzmxx9/FAVBEM+ePdtpY2+rhsfszI4dO0QAYlZWliiK6j7mpo73zJkzYmJionjw4EExJSVFfOutt5Tr1Hy8LWFmpg0sFgt2796NiRMnKpdpNBpMnDgRW7du9eDIOkZpaSkAICIiAgCwe/du1NTUOBx/nz59kJycrOrjnzVrFq655hqH4wJ883hXr16NESNGYMqUKYiJicHQoUPx8ccfK9dnZGQgLy/P4ZiNRiNGjRql2mMeM2YM1q5di2PHjgEA9u/fj82bN+Pqq68G4JvHLHPl2LZu3YqwsDCMGDFCuc3EiROh0Wiwffv2Th9zRygtLYUgCAgLCwPge8dst9tx5513Yu7cuejfv3+j633teOvrEhtNulthYSFsNhtiY2MdLo+NjcXRo0c9NKqOYbfbMXv2bIwdOxYDBgwAAOTl5cHf31/5QJDFxsYiLy/PA6Nsvy+++AJ79uzBzp07G13ni8d76tQpLFq0CHPmzMHTTz+NnTt34uGHH4a/vz+mTZumHJez17haj/nJJ5+EyWRCnz59oNVqYbPZ8Morr2Dq1KkA4JPHLHPl2PLy8hATE+NwvZ+fHyIiIlR//IBU9/bEE0/g9ttvVzZe9LVjfu211+Dn54eHH37Y6fW+drz1MZihZs2aNQsHDx7E5s2bPT2UDpOdnY1HHnkEv/zyCwwGg6eH0ynsdjtGjBiBV199FQAwdOhQHDx4EB988AGmTZvm4dF1jC+//BLLli3D8uXL0b9/f+zbtw+zZ89GQkKCzx4zSWpqanDrrbdCFEUsWrTI08PpELt378bbb7+NPXv2QBAETw+n03GaqQ2ioqKg1WobrWbJz89HXFych0blfg8++CC+//57rFu3Dt26dVMuj4uLg8ViQUlJicPt1Xr8u3fvRkFBAYYNGwY/Pz/4+flhw4YNeOedd+Dn54fY2FifOl4AiI+PR79+/Rwu69u3L06fPg0AynH50mt87ty5ePLJJ3Hbbbdh4MCBuPPOO/Hoo49i/vz5AHzzmGWuHFtcXFyjBQxWqxXnz59X9fHLgUxWVhZ++eUXJSsD+NYxb9q0CQUFBUhOTlY+x7KysvD//t//Q2pqKgDfOt6GGMy0gb+/P4YPH461a9cql9ntdqxduxajR4/24MjcQxRFPPjgg1i1ahV+++03pKWlOVw/fPhw6HQ6h+NPT0/H6dOnVXn8EyZMwIEDB7Bv3z7lNGLECEydOlU570vHCwBjx45ttNz+2LFjSElJAQCkpaUhLi7O4ZhNJhO2b9+u2mOurKyERuP4kafVamG32wH45jHLXDm20aNHo6SkBLt371Zu89tvv8Fut2PUqFGdPmZ3kAOZ48eP49dff0VkZKTD9b50zHfeeSf++OMPh8+xhIQEzJ07Fz/99BMA3zreRjxdgaxWX3zxhajX68UlS5aIhw8fFv/2t7+JYWFhYl5enqeH1m7333+/aDQaxfXr14u5ubnKqbKyUrnNfffdJyYnJ4u//fabuGvXLnH06NHi6NGjPThq96q/mkkUfe94d+zYIfr5+YmvvPKKePz4cXHZsmViYGCg+Nlnnym3+cc//iGGhYWJ3377rfjHH3+IkydPFtPS0sSqqioPjrztpk2bJiYmJorff/+9mJGRIX799ddiVFSU+Pjjjyu3UfMxl5WViXv37hX37t0rAhDffPNNce/evcrKHVeO7aqrrhKHDh0qbt++Xdy8ebPYq1cv8fbbb/fUIbWouWO2WCzi9ddfL3br1k3ct2+fw2eZ2WxWHkNNx9zSv3FDDVcziaK6jrc1GMy0w7vvvismJyeL/v7+4oUXXihu27bN00NyCwBOT4sXL1ZuU1VVJT7wwANieHi4GBgYKN54441ibm6u5wbtZg2DGV883u+++04cMGCAqNfrxT59+ogfffSRw/V2u1187rnnxNjYWFGv14sTJkwQ09PTPTTa9jOZTOIjjzwiJicniwaDQezevbv4zDPPOHyxqfmY161b5/R9O23aNFEUXTu2oqIi8fbbbxeDg4PF0NBQccaMGWJZWZkHjsY1zR1zRkZGk59l69atUx5DTcfc0r9xQ86CGTUdb2sIoliv/SURERGRyrBmhoiIiFSNwQwRERGpGoMZIiIiUjUGM0RERKRqDGaIiIhI1RjMEBERkaoxmCEiIiJVYzBDREREqsZghoi8xrlz53D//fcjOTkZer0ecXFxmDRpEn7//XcAgCAI+Oabbzw7SCLyOn6eHgARkezmm2+GxWLB0qVL0b17d+Tn52Pt2rUoKiry9NCIyItxOwMi8golJSUIDw/H+vXrMX78+EbXp6amIisrS/k7JSUFmZmZAIBvv/0WL774Ig4fPoyEhARMmzYNzzzzDPz8pN9rgiDg/fffx+rVq7F+/XrEx8djwYIFuOWWWzrl2IioY3GaiYi8QnBwMIKDg/HNN9/AbDY3un7nzp0AgMWLFyM3N1f5e9OmTbjrrrvwyCOP4PDhw/jwww+xZMkSvPLKKw73f+6553DzzTdj//79mDp1Km677TYcOXKk4w+MiDocMzNE5DX++9//4t5770VVVRWGDRuG8ePH47bbbsOgQYMASBmWVatW4YYbblDuM3HiREyYMAFPPfWUctlnn32Gxx9/HDk5Ocr97rvvPixatEi5zUUXXYRhw4bh/fff75yDI6IOw8wMEXmNm2++GTk5OVi9ejWuuuoqrF+/HsOGDcOSJUuavM/+/fvx97//XcnsBAcH495770Vubi4qKyuV240ePdrhfqNHj2ZmhshHsACYiLyKwWDAFVdcgSuuuALPPfcc/vrXv2LevHmYPn2609uXl5fjxRdfxE033eT0sYjI9zEzQ0RerV+/fqioqAAA6HQ62Gw2h+uHDRuG9PR09OzZs9FJo6n7iNu2bZvD/bZt24a+fft2/AEQUYdjZoaIvEJRURGmTJmCu+++G4MGDUJISAh27dqFBQsWYPLkyQCkFU1r167F2LFjodfrER4ejueffx7XXnstkpOTccstt0Cj0WD//v04ePAgXn75ZeXxV65ciREjRmDcuHFYtmwZduzYgX//+9+eOlwiciMWABORVzCbzXjhhRfw888/4+TJk6ipqUFSUhKmTJmCp59+GgEBAfjuu+8wZ84cZGZmIjExUVma/dNPP+Hvf/879u7dC51Ohz59+uCvf/0r7r33XgBSAfB7772Hb775Bhs3bkR8fDxee+013HrrrR48YiJyFwYzROTznK2CIiLfwZoZIiIiUjUGM0RERKRqLAAmIp/H2XQi38bMDBEREakagxkiIiJSNQYzREREpGoMZoiIiEjVGMwQERGRqjGYISIiIlVjMENERESqxmCGiIiIVI3BDBEREana/wdf021AyCBL+gAAAABJRU5ErkJggg==" }, "metadata": {}, "output_type": "display_data" }, { "name": "stdout", "output_type": "stream", "text": [ "Simulated mean: 2.004, True mean: 2\n", "Simulated error: 0.567, True error: 0.577\n" ] }, { "data": { "text/plain": [ "
" ], "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjcAAAHKCAYAAAD/zGr0AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy80BEi2AAAACXBIWXMAAA9hAAAPYQGoP6dpAABaw0lEQVR4nO3dd3xT9f7H8VeSLmaZLatQQNmjMi2KiKBVEERFcV0QFfUyRNHrBRf6cwCuiwNFUXAiKEMUEFBkiEwphbJlI9BSVgsF2iY5vz8CldoWkjbtSZP38/HIIyen5+S8Q0L7yff7PedrMQzDQERERMRPWM0OICIiIuJNKm5ERETEr6i4EREREb+i4kZERET8ioobERER8SsqbkRERMSvqLgRERERv6LiRkRERPyKihsRERHxKypuRERExK+ouBERERG/ouJGitxnn32GxWJhz549AXXsi3nxxRexWCxe2784X2dexzqf58iRI6Yc3xetWbOGDh06UKZMGSwWCwkJCWZH8plMxfl58YaS8pmTv6m4kUJJTEykd+/e1KlTh7CwMGrWrMn111/Pe++9Z3a0Qlu+fDkvvvgiJ06cMDtKkfHl1+jL2S4lKyuLO+64g2PHjvG///2PL7/8kjp16vhkppL87yySHxU3UmDLly+nTZs2rF+/ngEDBvD+++/z0EMPYbVaeeedd7K3+9e//sWZM2dM/+XuqeXLl/PSSy+ViF/6Bf03LshrLK73M79sJeHztHPnTvbu3ctTTz3Fww8/zH333UfFihV9MlNJ+pybpSR85iSnILMDSMn16quvEh4ezpo1a6hQoUKOnx0+fDh72WazYbPZijldYCmOf+P09HTKlClj+vtp9vHdcf7z/8//F2Yq7kznPy9m8laGkvCZk5zUciMFtnPnTpo2bZrnL8uIiIjs5YuN0di+fTv33Xcf4eHhVK1aleeffx7DMNi/fz+33HIL5cuXp1q1arz11ls5nv/+++8nOjo613HdGcuyd+9eBg4cSMOGDSlVqhSVK1fmjjvuyJXvP//5DwB169bFYrHkeg0HDhzggQceIDIyktDQUJo2bcrEiRNzHW/ZsmW0bduWsLAw6tevz0cffXTRfAXZP69/45MnT/L4448THR1NaGgoERERXH/99cTHx1/yNZ7/d9y8eTP33HMPFStW5Oqrr873WOcdOXKEO++8k/Lly1O5cmWGDh3K2bNns3/u7vt2sWz5HX/dunXcdNNNlC9fnrJly9KlSxdWrlyZ53F27NjB/fffT4UKFQgPD6d///6cPn06/zfBg+Pcf//9dOrUCYA77rgDi8XCtddee9HnrFevHvfdd1+u9Z07d85+rvy483nOL5O3PucX+7xcjLc+L5fKUNj3/WK/wwrzWZKio5YbKbA6deqwYsUKNm7cSLNmzQr0HH369KFx48aMHj2aOXPm8Morr1CpUiU++ugjrrvuOsaMGcPXX3/NU089Rdu2bbnmmmsKnXvNmjUsX76cu+66i1q1arFnzx4+/PBDrr32WjZv3kzp0qW57bbb2L59O9988w3/+9//qFKlCgBVq1YFIDk5mSuvvBKLxcLgwYOpWrUqP/30Ew8++CBpaWk8/vjjgGtM0g033EDVqlV58cUXsdvtjBw5ksjISLeyFmb/Rx99lGnTpjF48GCaNGnC0aNHWbZsGVu2bKFVq1aXfI3g+kN4+eWX89prr2EYxiWPeeeddxIdHc2oUaNYuXIl7777LsePH+eLL75w6/We5062C23atImOHTtSvnx5nn76aYKDg/noo4+49tprWbJkCe3bt8+Vs27duowaNYr4+Hg++eQTIiIiGDNmzEVzuXOcRx55hJo1a/Laa6/x2GOP0bZt24u+X6dOnWLPnj38+9//zvWzDRs2cM8991w0kzuf5/wyRUZGeuVzfp5Znxd3MxT0fb/Ua/D2c4oXGCIFtGDBAsNmsxk2m82IjY01nn76aWP+/PlGZmZmju0mTZpkAMbu3buz140cOdIAjIcffjh7nd1uN2rVqmVYLBZj9OjR2euPHz9ulCpVyujXr1/2un79+hl16tTJlen8817s2KdPn86134oVKwzA+OKLL7LXvfHGG7n2Pe/BBx80qlevbhw5ciTH+rvuussIDw/PPkavXr2MsLAwY+/evdnbbN682bDZbIY7//3c3T+v1xkeHm4MGjToos+f32s8/+94991359rnYu9nz549c2w7cOBAAzDWr19vGIb779vFsuV1/F69ehkhISHGzp07s9cdPHjQKFeunHHNNdfkOs4DDzyQ4zlvvfVWo3Llyrly/ZO7x1m0aJEBGN99990ln/P8Z2/+/Pk51u/fv98AjI8//vii+7v7ec4vkzc+5xf7vOSlKD4vF8tQ2Pf9Yp/5gj6nFC11S0mBXX/99axYsYKePXuyfv16Xn/9deLi4qhZsyY//PCDW8/x0EMPZS/bbDbatGmDYRg8+OCD2esrVKhAw4YN2bVrl1dylypVKns5KyuLo0ePctlll1GhQoXsLpuLMQyD6dOn06NHDwzD4MiRI9m3uLg4UlNTiY+Px+FwMH/+fHr16kXt2rWz92/cuDFxcXGXPE5h969QoQKrVq3i4MGDl9w2P48++qhH2w8aNCjH4yFDhgAwd+7cAme4FIfDwYIFC+jVqxf16tXLXl+9enXuueceli1bRlpaWo59/vm6OnbsyNGjR3NtV9jjuGPjxo0AtGzZMsf69evXA9CiRYuL7l/Yz3N+3P2cX8gXPi8Xy1CQ993T43njOaXwVNxIobRt25YZM2Zw/PhxVq9ezYgRIzh58iS9e/dm8+bNl9z/wj/aAOHh4YSFhWU3j1+4/vjx417JfObMGV544QWioqIIDQ2lSpUqVK1alRMnTpCamnrJ/VNSUjhx4gQff/wxVatWzXHr378/4Bq8mZKSwpkzZ7j88stzPUfDhg3dOk5h9n/99dfZuHEjUVFRtGvXjhdffNHjArFu3boebf/PrPXr18dqtRbp9UFSUlI4ffp0nv8mjRs3xul0sn///hzr//m5O38m08U+YwU5jjsSExOzu4gutGHDBqxW6yW7fAv7ec6Pu5/zC/nC5+ViGS72vmdmZpKUlJTj5nA4Lnm8gnyWpOhpzI14RUhICG3btqVt27Y0aNCA/v3789133zFy5MiL7pfXGQj5nZVgXNB/nt+gYXd+GQ0ZMoRJkybx+OOPExsbS3h4OBaLhbvuugun03nJ/c9vc99999GvX788t2nRooVbz1WU7rzzTjp27MjMmTNZsGABb7zxBmPGjGHGjBncdNNNbj3Hha0CBfHP96kw75s3ufMZKy4bN27M1WoDkJCQQL169S55tk9hP8/5cfdzfiFf+LxcLMPF3vfly5fTuXPnHOt3796d54Bmd59TzKPiRryuTZs2ABw6dKjIjlGxYsU8r8uxd+/eS+47bdo0+vXrl+MMrLNnz+Z6vvx+sVatWpVy5crhcDjo2rVrvsdxOByUKlWKP//8M9fPtm3bdsmcVatWLdT+4OoyGThwIAMHDuTw4cO0atWKV199Nbu4KcxVkvPy559/5vjmvGPHDpxOZ/YfCE/eN3ezVa1aldKlS+f5b7J161asVitRUVHuvQATjpOYmEifPn1yrHM6nfz6669uDaB39/Ocn8J+zgvDm5+XwmrZsiU///xzjnXVqlXz+nGkeKhbSgps0aJFeX47Od9f7k7XSUHVr1+f1NRUNmzYkL3u0KFDzJw585L72my2XLnfe++9XN8Gz39j/ucvV5vNxu2338706dOzx0tcKCUlJXu7uLg4vv/+e/bt25f98y1btjB//ny3chZ0f4fDkatLIiIigho1apCRkXHJ11hQ48aNy/H4/JWqzxdTnrxv7maz2WzccMMNzJo1K0d3RnJyMpMnT+bqq6+mfPnyBXk5RX6c892X//wi8O6773LkyBGaN2/uVi53Ps/5KeznvDC8+XkprIoVK9K1a9cct7CwMK889+nTp9m6dWuJmW7CH6jlRgpsyJAhnD59mltvvZVGjRqRmZnJ8uXLmTp1KtHR0dn98kXhrrvu4r///S+33norjz32GKdPn+bDDz+kQYMGlxxEefPNN/Pll18SHh5OkyZNWLFiBb/88guVK1fOsV3r1q0BePbZZ7nrrrsIDg6mR48elClThtGjR7No0SLat2/PgAEDaNKkCceOHSM+Pp5ffvmFY8eOAfDSSy8xb948OnbsyMCBA7Hb7bz33ns0bdo0xy/s/BR0/5MnT1KrVi169+5Ny5YtKVu2LL/88gtr1qzJ8Q0/v9dYULt376Znz57ceOONrFixgq+++op77rknu9vFk/fNk2yvvPIKP//8M1dffTUDBw4kKCiIjz76iIyMDF5//fUCv56iPk5iYiIACxYsYODAgTRq1IiVK1dmF69r165l1apVuU5lv5C7n+f8eONzXlDe/Lz4stWrV9O5c2dGjhzJiy++aHacwGDOSVriD3766SfjgQceMBo1amSULVvWCAkJMS677DJjyJAhRnJycvZ2FzuNMiUlJcdz9uvXzyhTpkyuY3Xq1Mlo2rRpjnULFiwwmjVrZoSEhBgNGzY0vvrqK7dOBT9+/LjRv39/o0qVKkbZsmWNuLg4Y+vWrUadOnVynG5uGIbx8ssvGzVr1jSsVmuu50lOTjYGDRpkREVFGcHBwUa1atWMLl265Dp1d8mSJUbr1q2NkJAQo169esb48ePzPJU1P+7s/8/XmZGRYfznP/8xWrZsaZQrV84oU6aM0bJlS+ODDz7I9fx5vcb83p/8/k3Pb79582ajd+/eRrly5YyKFSsagwcPNs6cOZNjf3fet4tly+v4hmEY8fHxRlxcnFG2bFmjdOnSRufOnY3ly5fn2Ca/15Xfc+bFneO4eyr4//73P8Nmsxlz5swx6tevb4SFhRnXX3+9kZiYaNSvX9+oVauWsXbt2os+h7uf54tlKuzn/GKfl7wUxeflYhkK+7578jssr23P/9uPHDnyoscR77EYhkY9iYiY4aGHHmLp0qVs377d7CgifkVjbkRETJKYmEiTJk3MjiHid1TciIiYwDAMNm/erOJGpAiouBERMcHu3bs5deqUihuRIqAxNyIiIuJX1HIjIiIifkXFjYiIiPiVgLuIn9Pp5ODBg5QrV87rl54XERGRomEYBidPnqRGjRpYrRdvmwm44ubgwYNemWdGREREit/+/fupVavWRbcJuOKmXLlygOsfxxvzzYiIiEjRS0tLIyoqKvvv+MUEXHFzviuqfPnyKm5ERERKGHeGlGhAsYiIiPgVFTciIiLiV1TciIiIiF9RcSMiIiJ+RcWNiIiI+BUVNyIiIuJXVNyIiIiIX1FxIyIiIn5FxY2IiIj4FRU3IiIi4ldU3IiIiIhfUXEjIiIifkXFjYiIiPgVFTciIiLiV1TciIiIiF8JMjuAiASu6OFzcjzeM7q7SUlExJ+o5UZERET8ioobERER8SsqbkRERMSvqLgRERERv6LiRkRERPyKihsRERHxKypuRERExK+ouBERERG/ouJGRERE/IqKGxEREfErKm5ERETEr6i4EREREb+i4kZERET8ioobERER8SsqbkRERMSvqLgRERERv6LiRkRERPyKihsRERHxKypuRERExK+ouBERERG/ouJGRERE/IqKGxEREfErKm5ERETEr6i4EREREb+i4kZERET8ioobERER8SsqbkRERMSvqLgRERERv6LiRkRERPyKihsRERHxKypuRERExK+ouBERERG/ouJGRERE/IqKGxEREfErKm5ERETEr6i4EREREb+i4kZERET8SpDZAUREilv08Dk5Hu8Z3d2kJCJSFNRyIyIiIn5FxY2IiIj4FdOLm3HjxhEdHU1YWBjt27dn9erVF91+7NixNGzYkFKlShEVFcUTTzzB2bNniymtiIiI+DpTi5upU6cybNgwRo4cSXx8PC1btiQuLo7Dhw/nuf3kyZMZPnw4I0eOZMuWLXz66adMnTqVZ555ppiTi4iIiK8ytbh5++23GTBgAP3796dJkyaMHz+e0qVLM3HixDy3X758OVdddRX33HMP0dHR3HDDDdx9992XbO0RERGRwGFacZOZmcnatWvp2rXr32GsVrp27cqKFSvy3KdDhw6sXbs2u5jZtWsXc+fOpVu3bsWSWURERHyfaaeCHzlyBIfDQWRkZI71kZGRbN26Nc997rnnHo4cOcLVV1+NYRjY7XYeffTRi3ZLZWRkkJGRkf04LS3NOy9AREREfJLpA4o9sXjxYl577TU++OAD4uPjmTFjBnPmzOHll1/Od59Ro0YRHh6efYuKiirGxCIiIlLcTGu5qVKlCjabjeTk5Bzrk5OTqVatWp77PP/88/zrX//ioYceAqB58+akp6fz8MMP8+yzz2K15q7VRowYwbBhw7Ifp6WlqcARERHxY6a13ISEhNC6dWsWLlyYvc7pdLJw4UJiY2Pz3Of06dO5ChibzQaAYRh57hMaGkr58uVz3ERERMR/mTr9wrBhw+jXrx9t2rShXbt2jB07lvT0dPr37w9A3759qVmzJqNGjQKgR48evP3221xxxRW0b9+eHTt28Pzzz9OjR4/sIkdEREQCm6nFTZ8+fUhJSeGFF14gKSmJmJgY5s2blz3IeN++fTlaap577jksFgvPPfccBw4coGrVqvTo0YNXX33VrJcgIiIiPsZi5Nef46fS0tIIDw8nNTVVXVQiJjNrAktNnClS8njy97tEnS0lIiIicikqbkRERMSvqLgRERERv6LiRkRERPyKihsRERHxKypuRERExK+ouBERERG/ouJGRERE/IqKGxEREfErKm5ERETEr6i4EREREb+i4kZERET8ioobERER8SsqbkRERMSvqLgRERERv6LiRkRERPyKihsRERHxKypuRERExK+ouBERERG/EmR2ABGRvEQPn5Pj8Z7R3U1KIiIljYobEZFCurAQUxEmYj51S4mIiIhfUXEjIiIifkXFjYiIiPgVFTciIiLiV1TciIiIiF/R2VIiUmz+eXq3iEhRUMuNiIiI+BW13IiIX9JFAEUCl4obESnZMk/DvuWQsh1S/wIMCC7FbdZ0Nhp12W7UAizqEhMJICpuRKTkMQzYvRRWfQQ7fgFHRq5N3g5x3e91RjDL2YHP7XEcJbyYg4qIGVTciEjJ8tda+OlpOPDH3+vCo6BmK6hQG6xBcDaNlauX09KykzrWwzxm/Z6HbXP4zBHHu/bbOE1YjqdUF5aIf1FxIyIlgz0Dfn4BVo13PQ4uDTH3Qpv+ENEELJYcm9+1bA6lOEsX6zoeCppLjHUnjwbN5mbbSp7OepjlzmYmvAgRKQ4qbkTE59WyHIZPukLSBteKlndD1xehXLWL7neGMGY7Y5mdeSXXWdfxf8GfUctyhK+CR/GW/Q4+cPTE0EmjIn5H/6tFxKc1s+xiZshIV2FTqhLc8y3cOv6ShU1OFn51tuL6jNeZYr8Wq8XgP8Hf8r/gDwjGXmTZRcQcKm5ExGe1sWxlasjLVLWkQmRzeHQZNIgr8POdIYzh9od5OmsAWYaNXrblTAh+i1AyvZhaRMym4kZEfFILy04mhbxBGUsGyxxNof9cCK/plef+1tGZAVlPcsYI4Vrbej4IfkctOCJ+RMWNiPie5M18ETKacpYzrHA04cGs/0BYea8eYrEzhn6Z/+WsEUwX2zreCv4QC06vHkNEzKHiRkR8S/oRmNyHCpZ04p2X8VDWk2QQUiSHWm005tGsJ8g0bPS0reDxoBlFchwRKV4qbkTEd9gz4du+kLqPPc5I+mc+TTqlivSQi50xjMgaAMDQoBl0s64s0uOJSNFTcSMivmPBs7D3dwgpx0NZT5JK2WI57HTnNUywdwPgreDxkLypWI4rIkVDxY2I+IQbrGtg9ceuB7d/wg6jVrEef7T9bpY6mlPKkgnTHoSsM8V6fBHxHhU3ImK6ahxlTPAE14MOj0HDG4s9gwMbj2cN4rBRAVK2wILnij2DiHiHihsRMZUFJ28Hf0hFyymoHgPXPW9almOU58msR10P1nwCW+ealkVECk7FjYiY6m7bIjrYNnPaCIXbP4Wgojkzyl2/OVtA7GDXg9lPwJkTpuYREc+puBER00RwnOFBkwF4w34nVLnM5ETnXPc8VL4MTiXBz+a1JIlIwai4ERHT/F/wZ5S3nCHBWZ/PHQWfVsHrgsOg53uu5fgvYNcSc/OIiEdU3IiIKeKsa7jRtoYsw8Z/swbg9LVfR3U6QJsHXcuzHwd7hqlxRMR9QWYHEJHAE0YGLwR/AcB4Rw+2GbUvuU/08Dk5Hu8Z3b1IsuXQ9UXYOhuO7eK1kUP52NGj+I4tIgXmY1+VRCQQPGybQ03LUf4yqvC+vZfZcfIXVh66jARgSND3VOWEuXlExC1quRGRYlWNozwa9CMAo7Lu8dq8Uf9s2fGalneTMPMtYqw7eTpoCv+xP1o0xxERr/G45WbSpEmcPn26KLKISAD4T/BUSlsyWONswBxne7PjXJrVyotZ/QC4I2gpLS07TA4kIpficXEzfPhwqlWrxoMPPsjy5cuLIpOI+KmWlh3cblsGwP9l9QUs5gZyU4JxGdMdHQEYEfwNGIbJiUTkYjwubg4cOMDnn3/OkSNHuPbaa2nUqBFjxowhKSmpKPKJiB95OmgqANMdHUk06pmcxjNvZN1JhhHMldYtsHOh2XFE5CI8Lm6CgoK49dZbmTVrFvv372fAgAF8/fXX1K5dm549ezJr1iycTmdRZBWRkmzXEq6ybSLTsPF2Vm+z03gsicp84bje9eCXl0C/50R8VqHOloqMjOTqq68mNjYWq9VKYmIi/fr1o379+ixevNhLEUWkxDMM+PVlACY7unCAqiYHKpgP7D05aZSCpA2w+Xuz44hIPgpU3CQnJ/Pmm2/StGlTrr32WtLS0pg9eza7d+/mwIED3HnnnfTr18/bWUWkpNo+H/5awxkjhHG+fOr3JRynPJ/Yu7keLHoVHHZzA4lInjwubnr06EFUVBSfffYZAwYM4MCBA3zzzTd07doVgDJlyvDkk0+yf/9+r4cVkZIlevgc6g7/kc1fPw3AZ444UqhgbqhC+sTRDUpXhqM7YMNUs+OISB48vs5NREQES5YsITY2Nt9tqlatyu7duwsVTET8ww3WtTSx7iXNKMV4ew+z4xRaOqWgw2Pwy0hY9ja0vMvsSCLyDx633HTq1IlWrVrlWp+ZmckXX7gup26xWKhTp07h04lICWcwMGgWAJ874kilrMl5vKTtgxBWwdV6s2mm2WlE5B88Lm769+9PampqrvUnT56kf//+XgklIv7hautGWlp3ccYIYZL9RrPjeE9oObhyoGv5t7ewoDOnRHyJx8WNYRhYLLkvvPXXX38RHh7ulVAi4h8G2VytNt84ruMY5U1O42XtH4HQ8nB4M9db15qdRkQu4PaYmyuuuAKLxYLFYqFLly4EBf29q8PhYPfu3dx4ox99MxORwtm3iljbZjINGxPs7s2iXWTzQ3nZ+ZxPBXVmcNAsBgd9z4LMNpSUKy6L+Du3i5tevXoBkJCQQFxcHGXL/t13HhISQnR0NLfffrvXA4pICfXbWwDMcHTkEJVNDlM0PrXfxAO2ebSw7ibWupkVzqZmRxIRPChuRo4cCUB0dDR9+vQhLCysyEKJSAmXtBH+nI/DsDDeUfLPkMrPccrzneMa+gX9zADbHBU3Ij7C4zE3/fr1U2EjIhe3YhwAPznbsceobnKYojXRcRNOw8J1tgQus/xldhwRwc3iplKlShw5cgSAihUrUqlSpXxvnho3bhzR0dGEhYXRvn17Vq9efdHtT5w4waBBg6hevTqhoaE0aNCAuXPnenxcESkiJ5Nh4zQAPnFzrE1JtteoxnxnGwAesul3kYgvcKtb6n//+x/lypXLXs7rbKmCmDp1KsOGDWP8+PG0b9+esWPHEhcXx7Zt24iIiMi1fWZmJtdffz0RERFMmzaNmjVrsnfvXipUqOCVPCLiBWs+AUcm1GpHwo7LzE5TLCbYu3OTbQ232pa5irtykWZHEglobhU3F84Tdf/993vt4G+//TYDBgzIvj7O+PHjmTNnDhMnTmT48OG5tp84cSLHjh1j+fLlBAcHA64xQCLiI7LOwB+fupZjB8IOc+MUl3ijAfHOy2hl3QFrJsB1z5kdSSSgeTzmJj4+nsTExOzHs2bNolevXjzzzDNkZma6/TyZmZmsXbs2e04qAKvVSteuXVmxYkWe+/zwww/ExsYyaNAgIiMjadasGa+99hoOhyPf42RkZJCWlpbjJiJFZMO3cPoohNeGRv47kDgvH9tvdi2s+QQyT5sbRiTAeVzcPPLII2zfvh2AXbt20adPH0qXLs13333H008/7fbzHDlyBIfDQWRkzubbyMhIkpKS8txn165dTJs2DYfDwdy5c3n++ed56623eOWVV/I9zqhRowgPD8++RUVFuZ1RRDxgGLDyQ9dy+4fB5vHUdSXaAmcb9jmrwpnj2WOORMQcHhc327dvJyYmBoDvvvuOTp06MXnyZD777DOmT5/u7Xw5OJ1OIiIi+Pjjj2ndujV9+vTh2WefZfz48fnuM2LECFJTU7Nvmq1cpIjsWgQpWyCkLLTqa3aaYufEypeO610PVn/sKvZExBQef7UyDAOn0zWPyi+//MLNN7uaYqOiorLPqHJHlSpVsNlsJCcn51ifnJxMtWrV8tynevXqBAcHY7PZstc1btyYpKQkMjMzCQkJybVPaGgooaGhbucSkQJaPQGASaev4qUXl5kcxhzfOq7l2VIzICkR9q+G2u3NjiQSkDxuuWnTpg2vvPIKX375JUuWLKF7d9epnrt3787VxXQxISEhtG7dmoULF2avczqdLFy4kNjY2Dz3ueqqq9ixY0d2cQWulqTq1avnWdiISDE5sR+2zwP4u/UiAKVSFpr3dj1Y/bG5YUQCmMfFzdixY4mPj2fw4ME8++yzXHaZ61TPadOm0aFDB4+ea9iwYUyYMIHPP/+cLVu28O9//5v09PTss6f69u3LiBEjsrf/97//zbFjxxg6dCjbt29nzpw5vPbaawwaNMjTlyEi3rT2MzCc/O5oyi6jhtlpzNV2gOt+8yw4ddjcLCIByuNuqRYtWuQ4W+q8N954I0d3kTv69OlDSkoKL7zwAklJScTExDBv3rzsFqB9+/Zhtf5df0VFRTF//nyeeOIJWrRoQc2aNRk6dCj//e9/PX0ZIuIt9kyI/xwI7FabbDVioFZb+GsNrP0cOv3H7EQiAafApzNkZmZy+PDhHF1EALVr1/boeQYPHszgwYPz/NnixYtzrYuNjWXlypUeHUNEvOefM3fvufcMpKdAuer8ktLKpFQ+pt3DruLmj4lw9RMBd+aYiNkKdLZUx44dKVWqFHXq1KFu3brUrVuX6Oho6tatWxQZRcSX/THRdd+qH/aCf1/yL01ugdJV4ORB2Dbn0tuLiFd5/Juof//+BAUFMXv2bKpXr+61qRhEpORpYNkPe38Hiw1a94P568yO5BuCQl3/Hr+95bqoX5NbzE4kElA8Lm4SEhJYu3YtjRo1Koo8IlKC3Gv7xbXQqBuUrwGouMnWuj/89jbsXgrHdkGlemYnEgkYHhc3TZo08eh6NiLin0pzltts565n0/ahYj/+P8f++JwKUXBZF9jxC8R/AV1fNDuRSMDweMzNmDFjePrpp1m8eDFHjx7VvE0iAaq7bSXlLGegUn2o28nsOL6p1blJhxMmgyPL3CwiAcTjlpvzE1126dIlx3rDMLBYLBedxFJE/Ecf22LXwhX3gZ+NvfNaq1CDG6FMVTiVDNvnQ+ObvfO8InJRHhc3ixYtKoocIlKC1LccoI11O3bDSlDMPWbH8V1BIRBzD/z+jqtrSsWNSLHwuLjp1EnNzyKB7k7bYgAWOWO4vlzec8HJOVf0dRU3O36G1AMQXtPsRCJ+z+MxNwC//fYb9913Hx06dODAgQMAfPnllyxbFpiT5YkEkiDs3Gb7DXBNFCmXUOUyqHM1GE5Y95XZaUQCgsfFzfTp04mLi6NUqVLEx8eTkZEBQGpqKq+99prXA4qIb+liXUdVSxopRjiLnDFmxykZWp8bWLzuS3BqXKJIUfO4uHnllVcYP348EyZMIDg4OHv9VVddRXx8vFfDiYjvOd8lNd1xja5I7K7GPSAsHFL3wy6NWxQpah4XN9u2beOaa67JtT48PJwTJ054I5OI+Kq0g1xrTQDgW4fG37ktuBS0uMu1vPZzc7OIBACPv3ZVq1aNHTt2EB0dnWP9smXLqFdPV+AU8WsJk7FZDFY7G7LLqGF2Gp+Ua2LR0d1dC636wuqPYNtPcPoYlK5kQjqRwOBxcTNgwACGDh3KxIkTsVgsHDx4kBUrVvDUU0/x/PPPF0VGEfGifP/4Xorz7wGxGkhcANWaQbXmkJQIG6dDuwFmJxLxWx4XN8OHD8fpdNKlSxdOnz7NNddcQ2hoKE899RRDhgwpiowi4gv2/g7Hd3PSKMUcR3uz05QYFxaTD9hieCE40XXFYhU3IkXG4zE3FouFZ599lmPHjrFx40ZWrlxJSkoKL7/8clHkExFfca7V5kfHlZwhzOQwJdMsRwewBsHBeDi81ew4In6rQNe5MQyDtLQ0IiMjadeuHWXLlvV2LhHxJRmnYMuPAEzTQOICO0o4XH6D68H6yeaGEfFjHhU3SUlJ9O3bl4oVKxIZGUlERAQVK1bkgQceIDk5uagyiojZts6GrHSoWJd443Kz05RsLe923W/4Vte8ESkibo+5SUtLo0OHDpw6dYr+/fvTqFEjDMNg8+bNfPPNNyxbtoz4+Hi14oj4o/XfuO5b3gWH/GuSzGLXIA5KVYSTh1zXvLmsq9mJRPyO28XNO++8g81mY9OmTVStWjXHz5577jmuuuoq3n33XZ555hmvhxQRE6UdhF1LXMst+sC8zebmKemCQqFZb1gzwTWwWMWNiNe53S01Z84cnnnmmVyFDUBERAQjRozgxx9/9Go4EfEBG74FDKgdC5Xqmp3GP8Sc65raOgfOppqbRcQPuV3cbN++nQ4dOuT78w4dOrBt2zavhBIRH2EYsH6Ka7lFH3Oz+JMaraBKQ7CfhU0zzU4j4nfcLm7S0tKoUKFCvj+vUKECaWlp3sgkIr4iaQOkbAFbKDTtZXYa/2GxQMw9ruWEb8zNIuKH3C5uDMPAas1/c4vFgmEYXgklIj5i/VTXfcMbXYNgxXta9AGLFfavhKM7s1dHD5+TfRORgnF7QLFhGDRo0ACLJe8zJVTYiPgZhx0Sv3Mtnz99WbynfHWo1xl2LnSNa+o8wuxEIn7D7eJm0qRJRZlDRHzNrkWQfhhKV9YZPUWlxZ2u4ibxW7h2uKu7SkQKze3ipl+/fkWZQ0R8zflr2zTrDbZgc7P4q0bdIagUHNsFB+KhVmuzE4n4hQJNvyAifu5smus0ZYCWOkuqyISWcxU44Gq9ERGvUHEjIrltnuU6TblKA9dpy1J0Wtzput843TXOSUQKze1uKRHxT/88K2fP6O6w4dxZUi36aBxIUat/nWtcU3oK7FpsdhoRv6CWGxHJKfUv2POba/l8q4IUHVswNL3VtayuKRGv8Li4WbRoUVHkEBFfsXG6677OVVChtrlZAkXzc0XkltmU4qy5WUT8gMfdUjfeeCO1atWif//+9OvXj6ioqKLIJSJmSZzmum92u7k5AklUO6hQB07s5XprPD84XVPd5NllKCKX5HHLzYEDBxg8eDDTpk2jXr16xMXF8e2335KZmVkU+USkGNW3HHBNuWANgia9zI4TOCyW7C7AW2y/mxxGpOTzuLipUqUKTzzxBAkJCaxatYoGDRowcOBAatSowWOPPcb69euLIqeIFIOethWuhfrXQZnK5oYJNOe6pq6xbqAimqdPpDAKNaC4VatWjBgxgsGDB3Pq1CkmTpxI69at6dixI5s2bfJWRhEpFgY9rMtdi816mxslEFVtANVbEmxx0N22yuw0IiVagYqbrKwspk2bRrdu3ahTpw7z58/n/fffJzk5mR07dlCnTh3uuOMOb2cVkSLUzLKbetYk1xVzG3XLXn/hRI6azLGInWu96aWuKZFC8bi4GTJkCNWrV+eRRx6hQYMGrFu3jhUrVvDQQw9RpkwZoqOjefPNN9m6dWtR5BWRInKL7VyrTcMbXVfOleLX7HachoU21u1EWZLNTiNSYnlc3GzevJn33nuPgwcPMnbsWJo1a5ZrmypVquiUcZESxIKTm20rXQ/UJWWe8tVZ7mwCwC3nuwhFxGMeFzcjR47kjjvuIDQ0NMd6u93O0qVLAQgKCqJTp07eSSgiRa6dZRvVLcdIM0rD5debHSegzXJeBZzvmjLMDSNSQnl8nZvOnTtz6NAhIiIicqxPTU2lc+fOOBwOr4UTkeJx/vTjnxzt6BMUeomtxZv+OY6pHO14OWgSl1kP0tSyh01GXZOSiZRcHrfcGIaBJY+5Zo4ePUqZMmW8EkpEik8wdm6yrQbgB2esyWnkJKX5xemarLSnTV1TIgXhdsvNbbfdBoDFYuH+++/P0S3lcDjYsGEDHTp08H5CESlSV1sTqWg5RYoRzgpnU7PjCPCDowM321bR07aC0fa7MTQNoIhH3C5uwsPDAVfLTbly5ShVqlT2z0JCQrjyyisZMGCA9xOKSJE63yU123ElTqw63dsHLHbGkGaUprrlGO0s21hlNDY7kkiJ4nZxM2nSJACio6N56qmn1AUl4gfCyOB661oAZjmuMjmNnJdJMHMd7bgraDG32H5nlV3FjYgnCnS2lAobEf/Q1RpPGUsGe50RJBj1zY4jFzh/1lQ32yqCsZucRqRkcavlplWrVixcuJCKFStyxRVX5Dmg+Lz4+HivhRORonX+wn0/OmOB/P9fS/Fb5WxMslGBSMsJrrGuZ6GztdmRREoMt4qbW265JXsAca9evYoyj4gUk/KcopM1AVCXlC9yYmW2I5YHg37iFttyFTciHnCruBk5cmSeyyJSct1oW0OIxcEWZxR/GrXMjiN5mOXowINBP3G9dS2lOWt2HJESQ+cXigSo85f3/9GhSzj4qg1GPXY7IyllyeR66x9mxxEpMdxqualYseJFx9lc6NixY4UKJCJFryrHibVuBnThPt9m4QfnVQy1zvh7YlMRuSS3ipuxY8cWcQwRKU49bCuxWgzWOi/nLyPi0juIaX5wxDI0aAbXWDdA+hEoU8XsSCI+z63ipl+/fkWdQ0SK0fnL+s9Sl5TP22nUJNEZTXPrHtj8PbR9yOxIIj7PreImLS2N8uXLZy9fzPntRMQ31bYkE2PdicOwMNdxZaGfT1c0LnqzHFe5ipvEaSpuRNzg9pib8zOBV6hQIc/xN+cn1NSs4CK+ree5gcS/O5txhHCT04g7Zjuu5JmgyVj3rYAT+6BCbbMjifg0t4qbX3/9lUqVKgGwaNGiIg0kIkXIMLK7pH5wqkuqpEiiMqucjYm1bYaN0+HqJ8yOJOLT3CpuOnXqlOeyiJQwyZtoYD1AhhHMfEdbs9OIB2Y5O7iKm8RpKm5ELsHtiTMvdPz4cT799FO2bNkCQJMmTejfv392646I+KjE7wBY5IzhJKVNDiOe+MnRjtGhn0PyRkjeDJFNzI4k4rM8vojf0qVLiY6O5t133+X48eMcP36cd999l7p167J06dKiyCgi3mAYsHEGoLOkSqJUysLl17sebJxmbhgRH+dxcTNo0CD69OnD7t27mTFjBjNmzGDXrl3cddddDBo0qCgyiog37F8Nqfs4ZYTxq/MKs9NIQTTv7bpP/M5VrIpInjzultqxYwfTpk3DZrNlr7PZbAwbNowvvvjCq+FExIvOfduf72xDBiEmhwlchTp1vsFNEFzGdcbUX2sgqp33gon4EY9bblq1apU91uZCW7ZsoWXLll4JJSJe5rDDppmA5pIq0UJKQ+ObXcvnxk+JSG5utdxs2LAhe/mxxx5j6NCh7NixgyuvdF0AbOXKlYwbN47Ro0cXTUoRKZw9SyE9BUpVYtnZZmankcJofgdsmOoqVuNGga1A54WI+DWLYVy649ZqtWKxWLjUpiXhIn5paWmEh4eTmpqqqylL4Ph+ECR8BW0eIHpZV7PTSAHtGd0dHFnwVkM4fRTumwGXdTE7lkix8OTvt1sl/+7du70STERMYM+ALT+6lpv1hmUnTI0jhWQLhqa3wppPXNe8UXEjkotbY27q1Knj9q0gxo0bR3R0NGFhYbRv357Vq1e7td+UKVOwWCz06tWrQMcVCQh//gwZqVC+JtSONTuNeEPzO1z3W36ErDPmZhHxQQXurN28eTP79u0jMzMzx/qePXt69DxTp05l2LBhjB8/nvbt2zN27Fji4uLYtm0bERER+e63Z88ennrqKTp27Fig/CIB4/w1UZreClaPzyEQH3L+TCsLTn4LrUKtzCOwfT407WVuMBEf43Fxs2vXLm699VYSExNzjMM5P5mmp2Nu3n77bQYMGED//v0BGD9+PHPmzGHixIkMHz48z30cDgf33nsvL730Er/99hsnTpzw9GWIBIaMU7Btnmv5/DVSpMQzsPKDowMDg35wnTWl4kYkB4+/xg0dOpS6dety+PBhSpcuzaZNm1i6dClt2rRh8eLFHj1XZmYma9eupWvXvwc4Wq1WunbtyooVK/Ld7//+7/+IiIjgwQcf9DS+SECIHj6H6OFzGPrSq2A/A5XqQ/UYs2OJF2VfZfrPBXDmhKlZRHyNxy03K1as4Ndff6VKlSpYrVasVitXX301o0aN4rHHHmPdunVuP9eRI0dwOBxERkbmWB8ZGcnWrVvz3GfZsmV8+umnJCQkuHWMjIwMMjIysh+npaW5nU+kpOtpO/cloXlvONe6Kv5hm1EbIprA4c2usTet/mV2JBGf4XHLjcPhoFy5cgBUqVKFgwcPAq5Bx9u2bfNuun84efIk//rXv5gwYQJVqlRxa59Ro0YRHh6efYuKiirSjCK+ogInucZ67hpVzdQl5ZcunI5BRLJ53HLTrFkz1q9fT926dWnfvj2vv/46ISEhfPzxx9SrV8+j56pSpQo2m43k5OQc65OTk6lWrVqu7Xfu3MmePXvo0aNH9jqn0+l6IUFBbNu2jfr16+fYZ8SIEQwbNiz7cVpamgocCQg32VYTbHGwyVmHplUbmB1HikKz22Hh/8HupXAyCcrl/r0pEog8brl57rnnsguK//u//2P37t107NiRuXPn8u6773r0XCEhIbRu3ZqFCxdmr3M6nSxcuJDY2NynrDZq1IjExEQSEhKybz179qRz584kJCTkWbSEhoZSvnz5HDeRQNDT6uqS+kHTLfivitEQ1R74e8Z3ESlAy01cXFz28mWXXcbWrVs5duwYFStWzD5jyhPDhg2jX79+tGnThnbt2jF27FjS09Ozz57q27cvNWvWZNSoUYSFhdGsWc5Lx1eoUAEg13qRQBbJMdpbXXPA/eiIZYTJeaQINb8D9q9ydU3FDjQ7jYhPKNSkJPv37wcoVDdPnz59SElJ4YUXXiApKYmYmBjmzZuXPch43759WHVtDhGP3GxbidVisMbZgIO4Nz5NSqgmveCn/8LBeDi6EyrXv+QuIv7O46rBbrfz/PPPEx4eTnR0NNHR0YSHh/Pcc8+RlZVVoBCDBw9m7969ZGRksGrVKtq3b5/9s8WLF/PZZ5/lu+9nn33G999/X6DjivirHrblgLqkAkLZqlC/s2s5cZq5WUR8hMfFzZAhQ/j44495/fXXWbduHevWreP111/n008/5bHHHiuKjCLiiaM7ibHuwm5Ymetof+ntpeQ7Px1D4ndw6bmQRfyex91SkydPZsqUKdx0003Z61q0aEFUVBR33303H374oVcDioiHzg0sXe5sylHCTQ4jxaJRdwgKg6N/wqH1UCPG7EQipvK45SY0NJTo6Ohc6+vWrUtISIg3MolIQRlG9lxSPzjVJRUwQstBgxtdy7rmjYjnxc3gwYN5+eWXc1z1NyMjg1dffZXBgwd7NZyIeCh5E6RsJcMIYr6jrdlppDid75raOB2cns3xJ+Jv3OqWuu2223I8/uWXX6hVqxYtW7YEYP369WRmZtKlSxfvJxQR951rtVnsjOEkpU0OI8Xq8ushNBxOHoK9y6FuR7MTiZjGreImPDxnv/3tt9+e47Gu+CviAwzD9a0dnSUVKKKHz8nxeE9sT1j3patrSsWNBDC3iptJkyYVdQ4RKay/1sCJfRBSloVnrzA7jZih+R2u4mbzLOj2BgSF5vjxhcXQntHdizudSLEp8NXxUlJSWLZsGcuWLSMlJcWbmUSkIM4PJG3YjbOEXnxb8U/RV0PZanD2BOxYeMnNRfyVx6eCp6enM2TIEL744ovsOaZsNht9+/blvffeo3Rp9fOLFDtH1t9zC7W4E9Zk5vjxP7svxE9Zba7JNFeOcxW7jbqZnUjEFB633AwbNowlS5bw448/cuLECU6cOMGsWbNYsmQJTz75ZFFkFJFL2bkITh+B0lWgXmez04iZmvd23W/7CTJOmptFxCQeFzfTp0/n008/5aabbsqeZbtbt25MmDCBadN06W8RU2yY6rpv3htshZoyTkq6GldApfpgPwNb55qdRsQUHhc3p0+fzp7U8kIRERGcPn3aK6FExAMZJ2HruW6nFneam0XMZ7HknI5BJAB5XNzExsYycuRIzp49m73uzJkzvPTSS8TGxno1nIi4Ycts17f0ypdBjVZmpxFfcL5rauevkH7E3CwiJvC4/Xrs2LHceOONuS7iFxYWxvz5870eUEQu4XyXVIs+rm/tIlUuh+oxcCgBNn8PbR8yOZBI8fK4uGnevDl//vknX3/9NVu3bgXg7rvv5t5776VUqVJeDygiF3EyCXYvcS2f/7YuAq6uqUMJkDhNxY0EHI+Km6ysLBo1asTs2bMZMGBAUWUSEXdtnA6GE2q1g0r1zE4jvqTZbbDgOdi3wnVxxwq1zU4kUmw8GnMTHBycY6yNiJgsu0tKA4nlH8rXcF3UD7Kn5RAJFB4PKB40aBBjxozBbrcXRR4RcdfhrXBoPViDoOltl95eAk/2WVO6TIcEFo/H3KxZs4aFCxeyYMECmjdvTpkyZXL8fMaMGV4LJyIXkfit6/6y66FMZXOziG9q0hPmPAnJGyF5s9lpRIqNx8VNhQoVcs0KLiLFzOmEDeeuYaIuKclPqYpw+Q2wbQ5snAa0NjuRSLHwuLjRDOEiPmD/SkjdByHloOFNZqcRX9a8t6u4SfwOaAXocgHi/9wec+N0OhkzZgxXXXUVbdu2Zfjw4Zw5c6Yos4lIfjac65Jq0hOCdQkGuYgGN0JIWTixj1aWP81OI1Is3C5uXn31VZ555hnKli1LzZo1eeeddxg0aFBRZhORvNgzYNNM17K6pORSQkpDo5sB6GlbbnIYkeLhdrfUF198wQcffMAjjzwCwC+//EL37t355JNPsFo9PulKRArqzwVw9gSUqw7RHc1OIz4kevicHI/3jO7uWmh+B2yYws22lbxs/xcObCakEyk+blcl+/bto1u3btmPu3btisVi4eDBg0USTETykfCN6755b7Dqj5S4oV4nKF2FKpY0rrJuNDuNSJFzu7ix2+2EhYXlWBccHExWVpbXQ4lIPtKPwJ/n5nCLudfcLFJy2IKh6a0A3KKuKQkAbndLGYbB/fffT2hoaPa6s2fP8uijj+a41o2ucyNShBK/A6cdalwBEY3NTiM+7sJuqlaWmswIhTjrGp7hQRNTiRQ9t4ubfv365Vp33333eTWMiFxCwteu+3+02vxzrIXIP8Ubl7PfWZUoawo3WP8AbjU7kkiRcbu40fVtREx2aAMkJYItBJrpQpriKQsznVfxmPV7etuWmh1GpEjpNCeRkmL9uYHEDW+C0pXMzSIl0jRHJwA6WhMh9S+T04gUHRU3IiWBI+vvC/e1vMfcLFJi7TMiWeVshNVi/F0si/ghFTciJcGfP8PpI1AmAi7rYnYaKcG+O9d6Q8JkMAxzw4gUEY/nlhIRE5wfSNziTrAFawCxFNhcR3teCvqMMsd2wb4VUKeD2ZFEvE4tNyK+Lv0IbJ/nWo5Rl5QUzmnCmO2IdT1Y97W5YUSKiIobEV+XOM11bZvqMRDZ1Ow04ge+c1zjWtg0EzJOmRtGpAiouBHxdesnu+51RWLxkj+MhlCpPmSlw+ZZZscR8ToVNyK+LHkTHFoP1mDXXFIiXmH5u4szQV1T4n9U3Ij4svgvXPcNb9S1bcS7Wt4NFivs/R2O7jQ7jYhXqbgR8VVZZ2H9FNdyq9zTn4gUSnhNqNfZtZww2dwsIl6m4kbEV235Ec6egPAoqH+d2WnEH11xbhzX+m/A6TA3i4gX6To3Ir4q/nPX/RX3gdVmbhbxTw27Q1gFSDsAOxfB5V3z3fSf11baM7p7EYcTKTi13Ij4oqM7Yc9vgMVV3IgUheAwaNHHtRz/malRRLxJLTcivuhcq80iRwv6j1rPntG1TA4kfqv1/bD6I9j2E5xMgnLVgNwtNSIliVpuRHyNIyt7gOcUh8baSBGLbAJR7V0Xilz3ldlpRLxCxY2Ir9n2E6SnkGKEs9B5hdlpJBC0vt91H/85OJ2mRhHxBhU3Ir7mXJfUd45O2NVzLMWhSS8IDYcT+2DXIrPTiBSaihsRX3JiH+xYCMBUx7XmZpHAEVIaWt7lWl47ydwsIl6g4kbEl6z7CjAguiN7jWpmp5FAcr5r6vzAYpESTMWNiK9wXDCg8/wfGpHicuHAYs03JSWcihsRX7H9J9fF1EpXhkY3m51GAtH5onrt51jQwGIpuVTciPiKNZ+47lv1dV1cTaS4ZQ8s3svV1o1mpxEpMBU3Ir7gyJ+wazFggdb9zU4jgeqCgcX32BaaHEak4FTciPiCNZ+67hvEQcU65maRwNbGVVxfb11LdY6aHEakYHQRDRGzZaZnX5GYtgPy3ESXwhdvy3cizIjGEN2RoD2/cU/QQt6y32lCOpHCUXEjUkzy/WOS+B1kpELFaKiv6RbEB7QbAHt+427br7xv70UGIWYnEvGIihsRMxnG3wOJ2zwIVvUUizkuLL5tWPgttBI1LMfoZl3FTGdHE5OJeE6/SUXM9NcaSEqEoDC44j6z04gA4MDGV/auAPQLmm9yGhHPqbgRMdPqCa77ZrdD6UrmZhG5wBTHdWQYQcRYd9HSssPsOCIeUXEjYpaTybD5e9dy2wdNjSLyT8coz2xnLAD9ghaYnEbEMypuRMyy5hNwZEKtdlCztdlpRHL53H4DAN2tK6lCqslpRNyn4kbEBKFkwh/nrm0TO9DcMCL52GDUZ53zMkItdu6y/Wp2HBG3qbgRMUEv2+9w+iiE14ZGPcyOI5Kv8603/wr6mWDsJqcRcY9OBRcpIvlfeM/gAdtPrsX2D4NN/w3Fd81xXskIYzKRlhP0sC5nhvMasyOJXJJabkSK2dXWjTS0/gXBZeCKf5kdR+Sisgjic3scAAOC5gKGuYFE3KDiRqSYPWib61q44j4oVcHULCLu+NrRhXQjlMbWfVyl2cKlBFB7uEgxqm85QGfbepyGhWt/a8S+pX93XWVPxyDiY1Ipy7eOa+kfNJ+HbXP43dnc7EgiF+UTLTfjxo0jOjqasLAw2rdvz+rVq/PddsKECXTs2JGKFStSsWJFunbtetHtRXzJA7Z5APzibMU+I9LkNCLum+i4EYdhoZNtAw0s+3P8LHr4nBw3EbOZXtxMnTqVYcOGMXLkSOLj42nZsiVxcXEcPnw4z+0XL17M3XffzaJFi1ixYgVRUVHccMMNHDhwoJiTi3imEmncZvsNgE/t3UxOI+KZ/UYk85xtAXjofNeqiI8yvbh5++23GTBgAP3796dJkyaMHz+e0qVLM3HixDy3//rrrxk4cCAxMTE0atSITz75BKfTycKFC4s5uYhn+gXNp5Qlk/XOeqwyGpkdR8Rjn9hdXae32H6Hk0kmpxHJn6nFTWZmJmvXrqVr167Z66xWK127dmXFihVuPcfp06fJysqiUqW85+XJyMggLS0tx02kuJXhDPfbXBMQfmjvCVjMDSRSAOuMy1njbECoxQ6rPzY7jki+TC1ujhw5gsPhIDIy59iDyMhIkpLc+1bw3//+lxo1auQokC40atQowsPDs29RUVGFzi3iqbttvxJuOc1OZ3XmO9uYHUekwM633pxYOp6mw6dpjI34JNO7pQpj9OjRTJkyhZkzZxIWFpbnNiNGjCA1NTX7tn///jy3EykqIWTxUJBrjMJ4Rw+Mkv3fTgLcz87W7HRWp4Ilnftsv5gdRyRPpv6WrVKlCjabjeTk5Bzrk5OTqVat2kX3ffPNNxk9ejQLFiygRYsW+W4XGhpK+fLlc9xEilMv2zKqWY5zyKjELMdVZscRKRQnVj6w3wLAQ0FzXfOkifgYU4ubkJAQWrdunWMw8PnBwbGxsfnu9/rrr/Pyyy8zb9482rRRE7/4LitOHrHNBuAT+01kEmxyIpHCm+XswH5nVapaUuljW2R2HJFcTG8fHzZsGBMmTODzzz9ny5Yt/Pvf/yY9PZ3+/fsD0LdvX0aMGJG9/ZgxY3j++eeZOHEi0dHRJCUlkZSUxKlTp8x6CSL5irOuob71ECeMMkxxXGd2HBGvsBPER46bAXgkaLYm1BSfY3px06dPH958801eeOEFYmJiSEhIYN68edmDjPft28ehQ4eyt//www/JzMykd+/eVK9ePfv25ptvmvUSRPJkwcmQoO8B+NxxA+mUMjeQiBd95+jEYaMCNS1H6WVbZnYckRx8YvqFwYMHM3jw4Dx/tnjx4hyP9+zZU/SBRLzgButamlj3ctIoxST7jWbHEfGqDEKYYO/Gs8GTGWibxQxHRxzYzI4lAvhIcSPid5xOhgbNAOAzRxwnKHfJXXRKrZQ0Xzu6MjDoB+pak+luXckPTg2YF99gereUiF/aOju71eYTTbUgfuo0YUw81yo5NGgGNhwmJxJxUXEj4m1OJywZA8AkRxyplDU5kEjRmeS4keNGWepbD3GL9Xez44gA6pYSKZS8upLirKv5KGQjaUYpTZApfu8UpfnIfjPDg6fweNB0fsjsYHYkEbXciHiTBSePnxtrM8lxo1ptJCB87riBFKM8ta0p3GFbYnYcERU3It7U3bqKxtZ9pBmlmGi/yew4IsXiDGHZVy0eEjQTss6anEgCnYobES8Jws5TQd8CrskF1WojgWSyowuHjErUsByD+M/NjiMBTmNuRC7in2Nq9ozunu+2d9kWEW1NJsUI5xOHxtpIYMkghPftvXg1eCIsfROuuA9CypgdSwKUWm5EvKA0ZxkaNB2Ad+y3cZq8Z6kX8WffOq5ln7MqpB/mrZceI3r4HF2/SUyh4kbECx6yzaWqJY09zkimODqbHUfEFFkE8Ya9DwCPBv1IVU6YG0gCloobkUKqRBoPB7lm/n7Tfid29fZKAPvRGUuCsx5lLBk8fq41U6S4qbgRKaQhQTMpazlLojOaOc72ZscRMZmF17LuBaCPbRGXWf4yOY8EIhU3IoVQ33KA+2y/ADDafjeG/kuJsNpozAJHa4IsToYHfWN2HAlA+k0sUlCGwQtBXxJscfCL4wp+dzY3O5GIzxhtvxu7YaWrbR3sXmp2HAkwGhwgcgGPzuzYPp9Otg1kGjZesd9XdKFESqBdRg0mO7rQN+hnmP8MPLwErDazY0mAUMuNSEHYM2H+CAAmOrqxx6huciAR3zPWfjupRmlISoQ/JpodRwKIihuRglj1IRzbRYoRzvvnLjsvIjkdozxv2u90Pfj1ZUg/Ym4gCRgqbkQ8dTIZlrwBwBj7XZyitMmBRHzX146uUK05nE2FX0aaHUcChIobEU8teA4yT0KNVkx3dDQ7jYhPc2KFbm+5Hqz7CvavMTeQBAQVNyKe2PkrJH4LFit0f1Onfou4o3Z7iHFd+4a5T4LTYW4e8Xv6zSziplAyYfYw14N2D0PN1uYGEilJur4EoeFwaD2s+dTsNOLnVNyIuGlI0Ew4vhvK1YDOz5odR6RkKVsVur7gWl74EpzYZ24e8WsqbkTccLnlLx6xueaPotvrEFbe3EAiJVHrB6B2LGSegtlPgGGYnUj8lIobkUuw4eCN4I8ItjigYTdodLPZkURKJqsVer4HtlDY8Qts+NbsROKnVNyIXMIjttnEWHeSZpSGbm+CxWJ2JJGSq8rl0Olp1/K8/8KpFHPziF9ScSNyEQ0t+3g8aBoAL2X1hfCaJicS8QNXDYXI5nDmOMx+XN1T4nUqbkTyEYSdt4LHE2Jx8LOjFdOduqaNiFfYgqHXOLAGw9bZkPC12YnEz2jiTJF8DA76nmbWPRw3yvJM1kOAuqNEvKZ6S+j8DCx8iVPfP8lN32ay34gEYM/o7iaHk5JOLTcieWhr2coQ20wAXsi6nxQqmBtIxB9dNRRqd6Cs5SxvB3+IFafZicRPqOVG5B/COcXYkHHYLAbTHR350dnB7EgiJVr08Dk5Hme3zFhtcOt4To5tT1vrdv5t+4Fxjl7FH1D8jlpuRHIweD34Y2pajrLLWY3ns/qbHUjEv1Wsw8isfgAMC/qO9pYtJgcSf6DiRuQC99l+Ic72B5mGjSFZQzhNmNmRRPzeDGdHpjs6YrMYvBfyHpw6bHYkKeFU3Iic99cfPB/0JQCj7PewyahrciCRQGHhuaz+bHfWJMJyAqY/qMk1pVBU3IiA65vi1H8RarEzz9GWSY4bzU4kElDOEMa/sx7ntBEKu5fC4tFmR5ISTAOKJaDkObDRkQXf9oOTB9nhrMFTWY+Q32nf/9xfRLxnp1GTZ7IeZGzIB7D0dajWDJrcYnYsKYHUciOy4DnYtxxCyvFw1jBOUdrsRCIB63vn1dD+364HMx+FQxvMDSQlkoobCWxrPoFV413Lt33ELqOGuXlEhPpLOrDU0RyyTsOUezT/lHhMxY0ErGut62Duf1wPOj8LjXRVVBFf4MDG4Kwh7HRWh9T9MPVeyDpjdiwpQVTcSEBqatnDuOB3wXBCzH1wzX/MjiQiF0ijLAOynoTQcNi/CqY/BA672bGkhFBxIwGnliWFiSGvU8aSAfWuhR5jwaJ5o0R8zS6jBtw9GWyhrgk25z6pGcTFLSpuJKBEcJyvg18l0nKCbc5acOcXrhmKRcQ3RV8Nt08ALLD2M50iLm7RqeASONKP8lXIa9SxHmavM4J/ZY5gdVi42alEAo7Hl1Rocgt0fxPmPAlLRkNoWegwpGjCiV9QcSOB4cwJ+OpWGlgPcMioxL1Zz3KYimanEhF3tX0ITh+HRa/Agud4ec5WPnV0+3sSznPynaRTAoqKG/F/6Ufhy16QtIEjRnnuyxzBX0ZVs1OJiKc6/Yd3ft7M0KAZPB/8FU4sRA83O5T4Io25Ef92Mgk+6wZJG6BMVe7LfIadRk2zU4lIAf3Pfjvv2XsBMDL4Sx61/WBuIPFJKm7Ef53YB5NugpStUK463D+XrUZts1OJSKFYeMt+R3aBMzx4CiOCvgZ0FpX8TcWN+KeD6+CTrnBsF1SoDf1/gqoNzE4lIl5h4S37nbySdS8AjwTN4c3gj7ChmcTFRcWN+J9t82BSNziVDBFNoP88qFTX7FQi4mWfOLrzZOaj2A0rvW1L+Sx4DOU5ZXYs8QEaUCz+wzBc80TNf8Z15eF6neHOz0Gne4v4renOaziRVYb3gt+no20jMy0jeSjrKXYb1YGcZ09d6swpnWnlP1TciH/ITIcfh0LidwBMsV/Lc5v7Y39xmcnBRKSoLXS2pnfmSCaEvEV96yG+D3meoVmDWeyMMTuamETdUlLyHd3pGl+T+B1YbLycdR/D7QOwq3YXCRibjWhuyXiFP5wNCLec5rOQ1/lP0BSNwwlQKm6k5DIMWPcVfHQNHN4MZSLg/tl86ugGaK4okUBzhHDuyXyWz+3XAzAo6Ae+CXmFahw1OZkUN321FZ+X16XaK3CSUcGfcJNtDQCrnI0YcnQIhz88XtzxRMSHZBLMSHt/VjkbMyZ4Au2s21gQ+l9eyuoLRjdNkhsgVNyIT3B/IJ9Bd+sqRgZ/QYTlBJmGjbftd/Cx42acaogUkXPmOq9kY2Zd3g1+nxjrTt4KGQ+Td0OPsVC+htnxpIipuJESo5blMC8HTaKzbT0AO5w1GJo1iE1GwU/z9ngCPxEpVoX5P7rPiOT2zBcZYJvDE0HTCP1zPoy7Eq57Dto8ADb9CfRX+qorPq8UZ3nMNoOfQ56ms209GUYQb2f1plvmqEIVNiLi/xzYGO/oSffM16BGK8hIhZ/+4xqrt/s3s+NJEVHZKr7L6YB1X7E49AUiLScAWO5ownP2B9hlqFlZRNy3w6gFD/0CayfBr6/A4U3w+c3QuKerJadqQ7MjihepuBGfY8XJY888y+CgmTSwHiDSAvucVXndfheznVeiM6FEpECsNmj7EDS9DRa9Cn9MhC0/wNbZ0PxO6ljasdeoZnZK8QIVN+IzgrFzq+03/m37gbrWZACOG2V5z34rXzm6kkmwyQlFxC+UrgTd3+KGZZczLGgaN9rWwIYpLAz5lh+dsXxi784mI/qST+PJ1Y+leKm4EfOl/sWwoG+5y7aYiHPdT8eMsky038QXjhtIo4yp8UTEP203ong06wma23cxLOg7OtvWc6vtd261/c5yRxPYZoXLb3C1+EiJouJGikWuU71fvQF2LoL4z2HbXB4LcgKQZFRkgr0b3zi6cJowM6KKSIBJNOrRP+u/NLfv4sGgudxsXUkH22b45i4oVwNi7oaYe6FyfbOjiptU3EixseCktWU7PW0r4K0hcPrvq4YudzThK0dXFjjbaNoEETFFolGPx7MGM4a76Rc0n0fLLYeTB+G3t1y3OldB01uh0c1QvrrZceUiLIZhGGaHKE5paWmEh4eTmppK+fLlzY7j/zLTYddiJn81gc629VS3HPv7Z2WqQrPboXV/ot/eaV5GEZE8BGOnizWePrZFXGPdgM1y/s+lBaLa8cquy1jsbMkOoyZ7Rt9satZA4MnfbxU34l32TDiUAHt/d11DYs8ycGRk//ikUYr5zrb07vc41O2UfREtXUxPRHxZNY6ysmea6+yq/aty/OyQUYnqV9wE9a+D6KuhnM64Kgqe/P1W+78UnGFA6n44tN51278K9q8B+5mc21Wow2dHGrLIeQUrnY3JIITel3UxJ7OISAEkURk69IUOgyHtIGyZzdLZX9DOutXVIp3wtesGEB4Ftdq6blHtILIZBGsMYXHyieJm3LhxvPHGGyQlJdGyZUvee+892rVrl+/23333Hc8//zx79uzh8ssvZ8yYMXTr1q0YEwcYw4BTh+Hon3DkTzi6A5I3ugqaM7knqjxqlGONsxGrnY1Y4mzBzqQa6No0IuI3yteA9g/Td2ZNQsmkjXUbX1+bDrsWQfIm15e+1P2waYZre4sVKtWHyCYQ0QQiGkPVxlCxDgSFmvta/JTpxc3UqVMZNmwY48ePp3379owdO5a4uDi2bdtGREREru2XL1/O3XffzahRo7j55puZPHkyvXr1Ij4+nmbNmpnwCszl/oSTF+F0QPoRSPsLUg9A6l+uW9pfcGI/HN3pumR5XqzBENGYKX9VYoNRn1XORuw0Ll3MqBtKRPxBBiH87mwON5z73ZtxEg7Ew19r/r6dPur6cnj0T9g864K9La5CqWL037fwKCgXCWWrubq3SlUkesTcHMfUNXUuzfQxN+3bt6dt27a8//77ADidTqKiohgyZAjDhw/PtX2fPn1IT09n9uzZ2euuvPJKYmJiGD9+/CWP529jbs4XCUHYKcNZ1v/3Ssg85foP9s/b2ROuIub0EUg/6voPd/oInD4GXOpjYIEKtaHK5VClgetS5dVbur6FBIWqWBERv/fPosKti/gZBpxKdrXoHN4Mh7e4lo/8CVnplz6oLYS/7OVJMSpw1ChHKmW4vUMzCKsApSqcu68IYeEQUsZ1Cy4NwaVcyzb/ufhpiRlzk5mZydq1axkxYkT2OqvVSteuXVmxYkWe+6xYsYJhw4blWBcXF8f3339flFEv7WQybP0RnE4wHOC0u1pEDMcF6xx/3zvtYDhzrrtwW0eWayCuPRPsZ8Ge4bp3XPg4g8TQdELJJMTicOV4p4D5LVbXN4XwWhBeE8rXdH2DCK/pak6tVI/o5xdC0oU7HTx3ExGRC+XZqn7BWMPo4bOpTBq1LYeJshymtuUwdSzJVLMcI8JygoZl0l3d/o5MalmOUMty5O8nW7XM/SDWYAgpDcFlXAVPcClXwWMLcd2sQeeWL1hnC/77sTXIdRFDixUs5+6ttnPLlpw/y162ulqdGvco7D9jgZla3Bw5cgSHw0FkZGSO9ZGRkWzdujXPfZKSkvLcPikpKc/tMzIyyMj4+2yd1FRX90paWlphouf21yaYMezS2xWBs+duANhCIaQshJaFkLL8kZTFSSOM05TipFGKuzrFQOnKrsuPl65ybrkylKpIs/9bCH/l9eybzt1ERAJX7Se+y/dn//yb4sw4fcl9UwgmhZqspWaun218Ks71JfbUYe4Z+yNVLamEW04RTjpPdox0tcSfTc15yzoNmWdc9zjPPVMmnMkETnj2YgurRmuo2cmrT3n+39idDifTx9wUtVGjRvHSSy/lWh8VFWVCmuJwEjiS708fefPH4osiIhIgwscW3/P931vePVbRWAyPhRfJM588eZLw8Is/t6nFTZUqVbDZbCQnJ+dYn5ycTLVqeV8noFq1ah5tP2LEiBzdWE6nk2PHjlG5cmUslpJ5Bk9aWhpRUVHs37/fL8YNlWR6L3yL3g/foffCt/jD+2EYBidPnqRGjRqX3NbU4iYkJITWrVuzcOFCevXqBbiKj4ULFzJ48OA894mNjWXhwoU8/vjj2et+/vlnYmNj89w+NDSU0NCcp9pVqFDBG/FNV758+RL7IfU3ei98i94P36H3wreU9PfjUi0255neLTVs2DD69etHmzZtaNeuHWPHjiU9PZ3+/fsD0LdvX2rWrMmoUaMAGDp0KJ06deKtt96ie/fuTJkyhT/++IOPP/7YzJchIiIiPsL04qZPnz6kpKTwwgsvkJSURExMDPPmzcseNLxv3z6sVmv29h06dGDy5Mk899xzPPPMM1x++eV8//33AXmNGxEREcnN9OIGYPDgwfl2Qy1evDjXujvuuIM77rijiFP5rtDQUEaOHJmru02Kn94L36L3w3fovfAtgfZ+mH4RPxERERFvsl56ExEREZGSQ8WNiIiI+BUVNyIiIuJXVNyIiIiIX1FxU8KMGzeO6OhowsLCaN++PatXrzY7UsBaunQpPXr0oEaNGlgsFvMnbw1go0aNom3btpQrV46IiAh69erFtm3bzI4VkD788ENatGiRfbG42NhYfvrpJ7NjCTB69GgsFkuOi+D6KxU3JcjUqVMZNmwYI0eOJD4+npYtWxIXF8fhw4fNjhaQ0tPTadmyJePGjTM7SsBbsmQJgwYNYuXKlfz8889kZWVxww03kJ6ebna0gFOrVi1Gjx7N2rVr+eOPP7juuuu45ZZb2LRJk++aac2aNXz00Ue0aNHC7CjFQqeClyDt27enbdu2vP/++4BrqoqoqCiGDBnC8OHDTU4X2CwWCzNnzsyeRkTMlZKSQkREBEuWLOGaa64xO07Aq1SpEm+88QYPPvig2VEC0qlTp2jVqhUffPABr7zyCjExMYwdO9bsWEVKLTclRGZmJmvXrqVr167Z66xWK127dmXFihUmJhPxPampqYDrj6qYx+FwMGXKFNLT0/Od/0+K3qBBg+jevXuOvx/+zieuUCyXduTIERwOR/a0FOdFRkaydetWk1KJ+B6n08njjz/OVVddpWlZTJKYmEhsbCxnz56lbNmyzJw5kyZNmpgdKyBNmTKF+Ph41qxZY3aUYqXiRkT8yqBBg9i4cSPLli0zO0rAatiwIQkJCaSmpjJt2jT69evHkiVLVOAUs/379zN06FB+/vlnwsLCzI5TrFTclBBVqlTBZrORnJycY31ycjLVqlUzKZWIbxk8eDCzZ89m6dKl1KpVy+w4ASskJITLLrsMgNatW7NmzRreeecdPvroI5OTBZa1a9dy+PBhWrVqlb3O4XCwdOlS3n//fTIyMrDZbCYmLDoac1NChISE0Lp1axYuXJi9zul0snDhQvVlS8AzDIPBgwczc+ZMfv31V+rWrWt2JLmA0+kkIyPD7BgBp0uXLiQmJpKQkJB9a9OmDffeey8JCQl+W9iAWm5KlGHDhtGvXz/atGlDu3btGDt2LOnp6fTv39/saAHp1KlT7NixI/vx7t27SUhIoFKlStSuXdvEZIFn0KBBTJ48mVmzZlGuXDmSkpIACA8Pp1SpUianCywjRozgpptuonbt2pw8eZLJkyezePFi5s+fb3a0gFOuXLlc487KlClD5cqV/X48moqbEqRPnz6kpKTwwgsvkJSURExMDPPmzcs1yFiKxx9//EHnzp2zHw8bNgyAfv368dlnn5mUKjB9+OGHAFx77bU51k+aNIn777+/+AMFsMOHD9O3b18OHTpEeHg4LVq0YP78+Vx//fVmR5MAouvciIiIiF/RmBsRERHxKypuRERExK+ouBERERG/ouJGRERE/IqKGxEREfErKm5ERETEr6i4EREREb+i4kZERET8ioobERER8SsqbkSkxFu2bBnBwcGcPXs2e92ePXuwWCzs3bvXxGQiYgYVNyJS4iUkJNC4cWPCwsKy161bt46KFStSp04dE5OJiBlU3IhIibd+/XquuOKKHOsSEhJo2bKlSYlExEwqbkSkxEtISCAmJibHunXr1uVaJyKBQcWNiJRoDoeDjRs35mq5iY+PV3EjEqBU3IhIibZt2zbOnj1LjRo1stetWLGCAwcOqLgRCVAqbkSkREtISADgvffe488//+Snn36ib9++AGRmZpqYTETMouJGREq0hIQE4uLi2LVrF82bN+fZZ5/lpZdeonz58rz77rtmxxMRE1gMwzDMDiEiUlBxcXG0bduWV155xewoIuIj1HIjIiXa+vXrad68udkxRMSHqLgRkRIrKSmJ5ORkFTcikoO6pURERMSvqOVGRERE/IqKGxEREfErKm5ERETEr6i4EREREb+i4kZERET8ioobERER8SsqbkRERMSvqLgRERERv6LiRkRERPyKihsRERHxKypuRERExK/8P7iJbJMVRt3pAAAAAElFTkSuQmCC" }, "metadata": {}, "output_type": "display_data" } ], "execution_count": 1 }, { "metadata": {}, "cell_type": "markdown", "source": "We want to compare the results of the MCMC with Least Squares (LS), see [Chapter 1](https://bayesian-statistics-for-astrophysics-2024.readthedocs.io/en/latest/lecture_notes/group1/Chapter%201.html). By implementing a simple LS program, we can find the best-fitting $\\,\\mu$ and the 68\\% confidence interval.", "id": "bf593cdef4965abd" }, { "metadata": { "ExecuteTime": { "end_time": "2024-12-19T18:26:35.562428Z", "start_time": "2024-12-19T18:26:33.282316Z" } }, "cell_type": "code", "source": [ "y = np.array([1, 2, 3])\n", "y_sigma = np.array([1, 1, 1])\n", "\n", "\n", "def chi_2(mu):\n", " return np.sum((y - mu[:, None]) ** 2 / (y_sigma ** 2), axis=1)\n", "\n", "\n", "mu_list = np.linspace(0, 4, 1000000)\n", "chi_2_list = chi_2(mu_list)\n", "min_chi_2 = min(chi_2_list)\n", "mu_fit = [x for x, y in zip(mu_list, chi_2_list) if y == min_chi_2][0]\n", "\n", "lower_bound = np.interp(-1 * (min_chi_2 + 1), -1 * chi_2_list[mu_list <= mu_fit], mu_list[mu_list <= mu_fit])\n", "upper_bound = np.interp(min_chi_2 + 1, chi_2_list[mu_list >= mu_fit], mu_list[mu_list >= mu_fit])\n", "\n", "plt.plot(mu_list, chi_2(mu_list), label=r'$\\chi^2(\\mu)$')\n", "plt.axhline(y=min_chi_2 + 1, color='red', linestyle='--', label='68% Confidence Level')\n", "plt.xlabel(r'$\\mu$')\n", "plt.ylabel(r'$\\chi^2$')\n", "plt.legend()\n", "plt.show()\n", "\n", "print(f'Calculated mean from LS: {mu_fit:.3f}±{(upper_bound - lower_bound) / 2:.3f}')" ], "id": "dcd0594d53b386dc", "outputs": [ { "data": { "text/plain": [ "
" ], "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjoAAAGxCAYAAABr1xxGAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy80BEi2AAAACXBIWXMAAA9hAAAPYQGoP6dpAABf/0lEQVR4nO3dd3iTVf8G8DtJ06QrKV20pYOyKbS0UIoMAWUpQ3EhyFJwvIoDecX5ivKKov4cuF63gIgDB1tlKXsU6KBljwKlg1IK6c48vz/SBgu0tKXtk3F/riuXkpy099OnSb495zznyIQQAkREREROSC51ACIiIqKmwkKHiIiInBYLHSIiInJaLHSIiIjIabHQISIiIqfFQoeIiIicFgsdIiIiclpuUgeQmsViQU5ODnx8fCCTyaSOQ0RERHUghEBxcTFCQ0Mhl9fcb+PyhU5OTg7Cw8OljkFEREQNkJWVhbCwsBofd/lCx8fHB4D1B6XRaCROQ0RERHVRVFSE8PBw2+d4TVy+0KkartJoNCx0iIiIHMy1pp1wMjIRERE5LRY6RERE5LRY6BAREZHTcvk5OkSNyWw2w2g0Sh2DyG4plUooFAqpY5ALYaFD1AiEEMjLy8PFixeljkJk93x9fREcHMy1y6hZsNAhagRVRU5QUBA8PT35Bk50FUIIlJWVIT8/HwAQEhIicSJyBSx0iK6T2Wy2FTn+/v5SxyGyax4eHgCA/Px8BAUFcRiLmhwnIxNdp6o5OZ6enhInIXIMVa8Vzmej5sBCh6iRcLiKqG74WqHmxEKHiIiInJZdFzqbN2/GqFGjEBoaCplMhmXLltXY9l//+hdkMhnmzZvXbPmIiIjIvtl1oVNaWopu3brhk08+qbXd0qVLsXPnToSGhjZTMiIiInIEdl3o3HrrrZgzZw7uuOOOGttkZ2fjiSeewOLFi6FUKpsxXe3KDWbsPVUodQyi65aVlYWBAwciOjoasbGx+Pnnn6WOREQOYu+pQpQZTJJmcOjLyy0WCyZOnIiZM2eiS5cudXqOXq+HXq+3/buoqKjRc2UVlmHYvM2wCIG9/xkCL5VD/5jJxbm5uWHevHmIi4tDXl4eevTogeHDh8PLy0vqaERkx8oMJoz/ahdkkGHt0/0R7ifNlal23aNzLW+99Rbc3Nzw5JNP1vk5c+fOhVartd3Cw8MbPVdYCw8E+qhQYbRg/cGzjf71iZpTSEgI4uLiAADBwcEICAhAYWH13srz588jKCgIJ0+erPPXHTt2LN59991GTEpE9mTDwXxUGC0I0qgQ1sJDshwOW+js3bsXH3zwARYsWFCvSxVfeOEF6HQ62y0rK6vRs8lkMoyMta74uWpfbqN/fSKp7N27F2az+Yo/EF5//XXcfvvtaN26dZ2/1n/+8x+8/vrr0Ol0jZySiOzB6srPvxExIZIuKeCwhc6WLVuQn5+PiIgIuLm5wc3NDadOncK///3vWt9sVSoVNBpNtVtTGNXNOjF60+FzKKrgoljk+AoLCzFp0iR88cUX1e4vKyvD119/jalTp9br63Xt2hVt27bFd99915gxicgOlOhN+PuwdauPEbHSbvXhsIXOxIkTsW/fPqSmptpuoaGhmDlzJtasWSN1PHRs6YN2Qd4wmC1Yt5/DV2S/7rnnHgQGBlYrYHbt2gV3d3esXbsWgHVu2+jRo/H888+jT58+1Z7/+++/Q6VS4YYbbrDdt3XrViiVSlRUVNjuO3nyJGQyGU6dOmW7b9SoUfjxxx+b6tCISCIbDp6F3mRBmwAvRIc0TYdCXdn1LNmSkhIcO3bM9u/MzEykpqbCz88PERERV+wrpFQqERwcjI4dOzZ31CtUDV/NW38UK/fl4K4eYVJHomYkhEC50SzJ9/ZQKurVTfzhhx9i9uzZ+O9//4uHH34YJSUlmDBhAh599FEMHToUQgjcf//9uPnmmzFx4sQrnr9lyxb06NGj2n2pqano3Lkz1Gq17b6UlBS0aNECkZGRtvsSExPx+uuvQ6/XQ6VSNeBoicgeVU3bGBEr7bAVYOeFzp49e3DTTTfZ/j1jxgwAwOTJk7FgwQKJUtXdyNhQzFt/FFuPFuBCqQEtvNyljkTNpNxoRvQsaXoWD/x3GDzd6/7SDgkJwfTp0/H555/j/PnzmDlzJlQqFd566y0AwLZt2/DTTz8hNjbWtmjnokWLEBMTAwA4derUFWtYpaWlIT4+vtp9qamp6NatW7X7QkNDYTAYkJeXV60AIiLHVVRhxKbD5wBIP2wF2HmhM3DgQAgh6ty+Pld8NId2Qd7oHKLBwdwirNmfh7GJEVJHIrqqDh06wNPTE7NmzcLixYuRlJRk643p168fLBZLjc8tLy+v1nMDWIua++67r9p9KSkptqu3qlTtZF1WVtYIR0FE9mD9gbMwmC1oF+SNji19pI5j34WOMxgZG4KDuUVYtS+XhY4L8VAqcOC/wyT73vUll8sRExOD//3vf3j77bev6HmpTUBAAC5cuGD7t9lsRkZGxhU9OsnJybjrrruq3Vd1mXpgYGC9MxORfaoathou8dVWVVjoNLFRsaH4vzWHsf14Ac4V6xHow3kIrkAmk9Vr+EhqVT2n3bt3x7///e96PTc+Pr7alVOHDx9GRUVFteGsHTt2IDs7+4oenYyMDISFhSEgIKDh4YnIblwoNWDzEeuw1W3dpB+2Ahz4qitHEeHviW5hWlgE8GcG19Qh+zRv3jzs2rULFosFcnn93haGDRuG/fv323p1UlNTAQAfffQRjh49ij/++AOTJk0CABgMhmrP3bJlC4YOHXr9B0BEduHP/XkwWQQ6h2jQLkj6YSuAhU6zGBlr/ct2JRcPJDuUnp6OF154AY899hgOHDgAk6l++9LExMSge/fuWLJkCQBroTNs2DCcOHECMTExeOmllzB79mxoNBp8+OGHtudVVFRg2bJleOihhxr1eIhIOivTcgAAo+ykNwdgodMsqmad7z5ZiDxdxTVaEzWfiooK3Hfffbj33nsxZ84cGAwGHDp0qN5fZ9asWfjggw9gsViQlpaGhIQErF69GhUVFUhOTsZ9990HnU6HRYsW2Z4zf/58JCYmVlt/h4gcV35RBXacOA/AOm3DXrDQaQahvh7oEdkCQgCr09mrQ/bj+eefR2lpKT7++GPbGjfz5s1DTk5Ovb7OiBEj8PDDDyM7OxtpaWm2S89ro1Qq8dFHHzU0OhHZmdXpuRACiI/wlWwDz6thodNMRtn2vqrfBwhRU1m7di0++eQTfPfdd/DxsY6l/+c//8GyZcswbdq0en+96dOnQ6lU4uzZs3UqdB588EG7WNyTiBqHbdjKjnpzAF511WyGx4Rg9qoDSDl9EVmFZXZV7ZJrGjp0KIzG6vuwPfjgg3jwwQcb/DWDg4PrtfYVETmHrMIyJJ++CJkMtk2t7QV7dJpJkEaNXlF+ADh8RUREzqXqc+2GKH8EadTXaN28WOg0o6odzTl8RUREzmRFatXVVvY1bAWw0GlWt3YNgUIuQ0Z2ETILSqWOQ0REdN2O5ZfgQG4R3OQy3No1WOo4V2Ch04z8vNzRp611x/VVaezVISIix1c1SnFj+wC73LyahU4zq5qNvoqLBxIRkYMTQmBFmv0OWwEsdJrdsC7BUCpkOHy2GEfOFksdh4iIqMEO5BbhxLlSqNzkGBLdUuo4V8VCp5lpPZXo3966U/NKDl8REZEDW5lmHZ24uVMQfNRKidNcHQsdCdwWZ+3eW56awzVHiIjIIQkh/rG3lX0OWwEsdCQxJLolPN0VOF1YhtSsi1LHISIiqrfk0xeRfbEcXu4K3NwpSOo4NWKhIwFPdzcMrRzLXJ7K4SsiInI8Vb05Q7sEQ61USJymZix0JHJ7XCsA1svyTGaLxGmIXM8XX3yB8PBwyOVyzJs3D6+++iri4uJqfc7999+P0aNHN0s+R7Zx40bIZDJcvHhR6ijURMwWYVsNeVQ3+9ry4XIsdCTSr30AWngqUVBiwPbj56WOQy4sOzsbEyZMgL+/Pzw8PBATE4M9e/bYHi8pKcHjjz+OsLAweHh4IDo6Gp999lm1rzFjxgz4+fkhPDwcixcvrvbYzz//jFGjRtUpi8FgwNtvv41u3brB09MTAQEB6Nu3L+bPn3/FvlzXo6ioCI8//jiee+45ZGdn4+GHH8YzzzyDDRs2NNr3kFLr1q0xb948qWOQE9uVeR7nivXQeijRr12g1HFqxU09JaJUyDEiNgTf7TyN5ak56N/Bvn9RyDlduHABffv2xU033YQ//vgDgYGBOHr0KFq0aGFrM2PGDPz111/47rvv0Lp1a6xduxaPPfYYQkNDcdttt2HlypX4/vvvsXbtWhw9ehRTpkzBsGHDEBAQAJ1Oh5deegnr16+/ZhaDwYBhw4YhLS0Nr732Gvr27QuNRoOdO3finXfeQXx8/DV7XOrq9OnTMBqNGDFiBEJCLv016u3t3Shfn8jZVQ1bDY8JhrubffeZ2Hc6Jze6cvhqzf48VBjNEqehJlFaWvOtoqLubcvL69a2nt566y2Eh4dj/vz5SExMRFRUFIYOHYq2bdva2mzfvh2TJ0/GwIED0bp1azz88MPo1q0bkpKSAAAHDx7EwIEDkZCQgHHjxkGj0SAzMxMA8Oyzz+LRRx9FRETENbPMmzcPmzdvxoYNGzBt2jTExcWhTZs2uO+++7Br1y60b98eAKDX6/Hkk08iKCgIarUa/fr1w+7du21fp2rYZMOGDUhISICnpyf69OmDw4cPAwAWLFiAmJgYAECbNm0gk8lw8uTJK4auzGYzZsyYAV9fX/j7++PZZ5+94ipJi8WCuXPnIioqCh4eHujWrRt++eWXOmepsnLlSvTs2RNqtRoBAQG44447bI/p9Xo888wzaNWqFby8vNCrVy9s3Ljxmj/P2ixfvhzdu3eHWq1GmzZtMHv2bJhMJgDAfffdh3vvvbdae6PRiICAAHz77bd1Om5ybgaTBb+n5wG4tAiuXRMuTqfTCQBCp9M1+/c2my2iz9wNIvK5VWJVWk6zf39qHOXl5eLAgQOivLz8ygeBmm/Dh1dv6+lZc9sBA6q3DQi4ert66ty5s5g+fbq4++67RWBgoIiLixNffPFFtTYPPfSQSEhIEGfOnBEWi0X89ddfwtvbW2zatEkIIcSff/4p2rZtKwoLC8WePXuEj4+PKCwsFFu2bBEJCQnCZDLVKUtsbKwYOnToNds9+eSTIjQ0VPz+++9i//79YvLkyaJFixbi/PnzQggh/v77bwFA9OrVS2zcuFHs379f3HjjjaJPnz5CCCHKysrE+vXrBQCRlJQkcnNzhclkEq+88oro1q2b7fu89dZbokWLFuLXX38VBw4cEFOnThU+Pj7i9ttvt7WZM2eO6NSpk/jzzz/F8ePHxfz584VKpRIbN26sUxYhhFi1apVQKBRi1qxZ4sCBAyI1NVW88cYbtscffPBB0adPH7F582Zx7Ngx8X//939CpVKJI0eO1PgzioyMFO+///5VH9u8ebPQaDRiwYIF4vjx42Lt2rWidevW4tVXX7Xl8fDwEMXFxbbnrFy5Unh4eIiioqJ6HfeFCxeumqHW1wzZvTUZuSLyuVUi8fV1wmS2SJajrp/fLHQkLHSEEOLNPw6KyOdWiYcW7pbk+9P1c+RCR6VSCZVKJV544QWRnJwsPv/8c6FWq8WCBQtsbSoqKsSkSZMEAOHm5ibc3d3FwoULq32dV155RbRt21Z07dpV/Pbbb0Kv14uuXbuKPXv2iI8++kh06NBB9OnTR2RkZNSYxcPDQzz55JO15i0pKRFKpVIsXrzYdp/BYBChoaHi7bffFkJc+pBdv369rc3q1asFANs5SklJEQBEZmZmtWP4Z6ETEhJi+5pCCGE0GkVYWJit0KmoqBCenp5i+/bt1TJOnTpVjBs3rs5ZevfuLcaPH3/V4z116pRQKBQiOzu72v2DBg0SL7zwQo0/p9oKnUGDBlUrpIQQYtGiRSIkJMR2nAEBAeLbb7+1PT5u3Dhx77331vu4Weg4p0e/2yMin1sl5qzaL2mOun5+c46OxG6PC8WnG49j4+Fz0JUZofW0z5UlqYFKSmp+THHZ5Zj5+TW3lV82ynzyZIMj/ZPFYkFCQgLeeOMNAEB8fDwyMjLw2WefYfLkyQCAjz76CDt37sSKFSsQGRmJzZs3Y9q0aQgNDcXgwYMBAK+++ipeffVV29edPXs2Bg8eDKVSiTlz5iA9PR2rVq3CpEmTsHfv3qtmEXVYPPP48eMwGo3o27ev7T6lUonExEQcPHiwWtvY2Fjb/1fNw8nPz6/TMJpOp0Nubi569eplu8/NzQ0JCQm2nMeOHUNZWRmGDBlS7bkGgwHx8fF1zpKamoqHHnroqjnS09NhNpvRoUOHavfr9Xr4+/tf8ziuJi0tDdu2bcPrr79uu89sNqOiogJlZWXw9PTEmDFjsHjxYkycOBGlpaVYvnw5fvzxx3ofNzmfogoj1h+0vldVXT1s71joSKxTsAadgn1wKK8Yf+7Pxb09r/0mTA7Ey0v6trUICQlBdHR0tfs6d+6MX3/9FQBQXl6OF198EUuXLsWIESMAWD+0U1NT8c4779gKnX86dOgQvvvuO6SkpOCbb75B//79ERgYiDFjxmDKlCkoLi6Gj4/PFc/r0KEDDh061CjHBVgLoCoymQyAtbBrLCWVRezq1avRqlX1N3yVSlXnLB4eHrV+D4VCgb1790JxWWHc0InTJSUlmD17Nu68884rHlOr1QCA8ePHY8CAAcjPz8e6devg4eGBW265xfZ8oG7HTc7nz4w8GEwWtA/yRpdQjdRx6oSTke1A1ZYQy1K4eCA1r759+14xMfbIkSOIjIwEYJ2EajQaIb+sR0mhUFy1aBBC4JFHHsF7770Hb29vmM1m22XhVf81m68+8f6+++7D+vXrkZKScsVjRqMRpaWlaNu2Ldzd3bFt27Zqj+3evfuKgu16aLVahISEYNeuXbb7TCZTtd6o6OhoqFQqnD59Gu3atat2Cw8Pr/P3io2NrfGy9vj4eJjNZuTn51/xPYKDgxt0bN27d8fhw4ev+Hrt2rWznec+ffogPDwcP/30ExYvXox77rnHVqw11nGTY1qemg0AGB3fyla02zv26NiBUbGhePvPw9iZeR55ugoEa9VSRyIX8fTTT6NPnz544403MGbMGCQlJeGLL77AF198AQDQaDQYMGAAZs6cCQ8PD0RGRmLTpk349ttv8d57713x9b766isEBgba1s3p27cvXn31VezcuRN//PEHoqOj4evre9Us06dPx+rVqzFo0CC89tpr6NevH3x8fLBnzx689dZb+PrrrxEXF4dHH30UM2fOhJ+fHyIiIvD222+jrKwMU6dObdSfzVNPPYU333wT7du3R6dOnfDee+9VWwDPx8cHzzzzDJ5++mlYLBb069cPOp0O27Ztg0ajsQ39Xcsrr7yCQYMGoW3bthg7dixMJhN+//13PPfcc+jQoQPGjx+PSZMm4d1330V8fDzOnTuHDRs2IDY21tbLdjXZ2dlITU2tdl9kZCRmzZqFkSNHIiIiAnfffTfkcjnS0tKQkZGBOXPm2Nred999+Oyzz3DkyBH8/fffjX7c5HjydBW2dd9us+O9ra7QHBOG7JnUk5Gr3PW/bSLyuVXiy83HJc1B9efoEytXrlwpunbtKlQqlejUqdMVV13l5uaK+++/X4SGhgq1Wi06duwo3n33XWGxVL/aIi8vT0RGRl4xcXb27NnCz89PdOrUSezatavWLBUVFWLu3LkiJiZGqNVq4efnJ/r27SsWLFggjEajEML6837iiSdEQECAUKlUom/fviIpKcn2Na42Efbyycd1mYxsNBrFU089JTQajfD19RUzZswQkyZNqnbVlcViEfPmzRMdO3YUSqVSBAYGimHDhtmuSKtLFiGE+PXXX0VcXJxwd3cXAQEB4s4777Q9ZjAYxKxZs0Tr1q2FUqkUISEh4o477hD79u2r8ecYGRkpAFxxW7RokRDCeqVcnz59hIeHh9BoNCIxMfGK837gwAEBQERGRl5xrhty3P/k6K8ZV/XFpuMi8rlV4q7/bZM6ihCi7p/fMiFce/vsoqIiaLVa6HQ6aDTSjTcu2nESLy/fj5hWWqx8op9kOaj+KioqkJmZiaioKNscByKqGV8zjmnEh1uwP6cIr43uiok3REodp86f35yjYydGxIbCTS5DerYOx8/VcqUOERFRMzuWX4z9OUVwk8swMsa+97a6HAsdO+Hn5Y4b2wcA4I7mRERkX6oulhnYMRAtvNwlTlM/LHTsSNWaBCtSs+u0pggREVFTE0JgWeXVVo6yds4/sdCxI0OiW8JDqcDJ82XYd0YndRwiIiLsPXUBZy6Uw8tdgcGdW0odp95Y6NgRL5UbhkRbf4mWpmRLnIbqi71wRHXD14pjqerNGdY1GB7uimu0tj8sdOzM6Hjr2gQr03JgNDfeKq7UdKoWUisrK5M4CZFjqHqt/HPFaLJPBpMFq/flAgDuiHe8YSuACwbanRvbByLA2x0FJQZsOXoON3dyvG5CV6NQKODr64v8yr2qPD09HWbFUKLmJIRAWVkZ8vPz4evre8W2FmR/Nh85hwtlRgT6qNCnbYDUcRqEhY6dUSrkGNUtFPO3ncRvydksdBxE1XL8+bVtzElEAABfX98Gb2FBzatq2GpUbCgUcsf8A46Fjh26Mz4M87edxNoDZ1FUYYRGze5deyeTyRASEoKgoCDbnk5EdCWlUsmeHAdRojdh/cGzAC5Nq3BELHTsUNdWGrQP8sbR/BL8kc4dzR2JQqHgmzgROYU1GXmoMFrQJtALMa20UsdpME5GtkMymQx3dLdO+votmVdfERFR86sathod5zg7lV8NCx07Zf3FAnZlFiKrkFfzEBFR88nTVWDbsQIAwO1xjjtsBbDQsVuhvh7o3cYfALA8lb06RETUfJalZsMigJ6tWyDS30vqONeFhY4dq1qz4LcUbglBRETNQwiBX/eeAQDc1T1M4jTXj4WOHbs1JgRqpRwnzpUijVtCEBFRM9ifU4Sj+SVwd5NjeKxj7VR+NXZd6GzevBmjRo1CaGgoZDIZli1bZnvMaDTiueeeQ0xMDLy8vBAaGopJkyYhJ8d5dv72VrlhWBfrWhNLk89InIaIiFzBL5W9OUOjWzrF8iZ2XeiUlpaiW7du+OSTT654rKysDMnJyXj55ZeRnJyM3377DYcPH8Ztt90mQdKmUzV8tXJfLgwmbglBRERNx2i2YEWatcPgrh6OP2wF2Pk6OrfeeituvfXWqz6m1Wqxbt26avd9/PHHSExMxOnTpxER4Rxrz/RrF4BAHxXOFeux6cg526afREREjW3T4XMoLDUgwFuFG9s55pYPl7PrHp360ul0kMlk8PX1rbGNXq9HUVFRtZs9c1PIcXs366V9S1M4fEVERE3n18ppEqPjQuGmcI4SwTmOAkBFRQWee+45jBs3DhqNpsZ2c+fOhVartd3Cw8ObMWXDVC0euP5gPnRl3F6AiIga38UyAzYctO7X5yzDVoCTFDpGoxFjxoyBEAKffvpprW1feOEF6HQ62y0rK6uZUjZcdIgGnYJ9YDBZsDo9V+o4RETkhFbty4XBbEHnEA06h9TcYeBoHL7QqSpyTp06hXXr1tXamwMAKpUKGo2m2s3eyWQy26RkDl8REVFTqBq2uqtyFMFZOHShU1XkHD16FOvXr4e/v7/UkZrM7ZVbQuw+eQGnz3NLCCIiajwnzpUg5fRFKOQy3ObgWz5czq4LnZKSEqSmpiI1NRUAkJmZidTUVJw+fRpGoxF333039uzZg8WLF8NsNiMvLw95eXkwGAzSBm8CwVo1+lXOgF+awi0hiIio8VR9rvRvH4AgH7XEaRqXXRc6e/bsQXx8POLj4wEAM2bMQHx8PGbNmoXs7GysWLECZ86cQVxcHEJCQmy37du3S5y8adxZ2Z34a/IZWCzcEoKIiK6fxSLwW7K10LnTCbZ8uJxdr6MzcODAWvd4crX9n27pEoKXVftxurAMSScLcUMb5x2qIyKi5rErsxDZF8vho3ZzyrXa7LpHh6rzcFdgZOW+I1VLdBMREV2P3yonIY+MDYFaqZA4TeNjoeNg7q5c2+D39FyU6k0SpyEiIkdWZjDh98plS5xhp/KrYaHjYHpEtkBUgBfKDGbbLycREVFDrN1/FqUGMyL8PNEjsoXUcZoECx0HI5PJbL06P3P4ioiIrkPV2jl3dm8FmUwmcZqmwULHAd3ZvRXkMiApsxCnzpdKHYeIiBxQ9sVybD1WAMB5h60AFjoOKUTrgX7tAwEAv7JXh4iIGuDXvWcgBNCnrT/C/TyljtNkWOg4qKrhq1+Ts7mmDhER1YvFIvDzXutej/ckOG9vDsBCx2ENjW4JH7Ubsi+WY8eJ81LHISIiB7IrsxBZheXwUbnhli4hUsdpUix0HJRaqcBt3az7kfy8x/53YCciIvtR9bkxslsoPNydb+2cf2Kh48DuSQgHAPy5Pw9FFUaJ0xARkSMoqjDi9wzr8iRjnHzYCmCh49C6hWnRLsgbFUYLVu/jmjpERHRtq/flosJoQbsgb8SF+0odp8mx0HFgMpkM91StqcPhKyIiqoMllZ8XYxLCnHbtnH9ioePg7ohvBYVchuTTF3H8XInUcYiIyI4dyy9GyumLUMhluCPe+YetABY6Di9Io8aADtY1dbjRJxER1ebnPdbPiZs6BiHQRyVxmubBQscJVA1f/ZZ8BmauqUNERFdhNFvwa3I2ANeYhFyFhY4TuLlzEHw9lThbpMeWo+ekjkNERHZo0+FzKCjRI8DbHTd1CpI6TrNhoeMEVG4KjI5rBeBStyQREdE/VU1CviO+FZQK1/n4d50jdXJjKtfUWXsgD+dL9BKnISIie1JQosdfh/IBXFqDzVWw0HES0aEaxIZpYTQLLE3JljoOERHZkWUp2TBZBLqF+6JDSx+p4zQrFjpOZGzPCADAj7uzIAQnJRMRESCEqLZ2jqthoeNERnULgYdSgWP5JUg+fUHqOEREZAf2ndHhyNkSqNzkGFW5R6IrYaHjRHzUSoyMte5C+2MSV0omIiLgp8renFu6BkOjVkqcpvmx0HEyYxOtk8xW7ctFMTf6JCJyaaV6E1ak5gAA7u3pWpOQq7DQcTLdI1qgXZA3yo1mrEzjRp9ERK5s9b5clOhNaO3vid5t/KWOIwkWOk5GJpNhbGXV/tPu0xKnISIiKf1Q+Tlwb88Il9jA82pY6Dgh62JQMqSd0eFATpHUcYiISAKH86wbeLrJZbi7h+tdbVWFhY4T8vdWYWh0MAD26hARuaofkqzv/4M7t3SZDTyvhoWOk6qadLY0JRsVRrPEaYiIqDlVGM22xWOrLlJxVSx0nFS/dgFo5euBogoT/szIkzoOERE1oz8z8qArN6KVrwdubB8odRxJsdBxUnK5zLb/1Y8cviIicilVw1ZjEsKhkLvmJOQqLHSc2D0JYZDJgJ0nCpFZUCp1HCIiagYnzpVgV2Yh5DJgTE/XnYRchYWOEwv19cCADtYuy6p9ToiIyLn9tNv6fj+wYxBCtB4Sp5EeCx0nV7Wmzi97z8BotkichoiImpLBZMEve88AuPT+7+pY6Di5mzu1RIC3O84V6/HXoXyp4xARURNad+AszpcaEOSjws2dgqSOYxdY6Dg5dzc57qpcKOr7XZyUTETkzKouPrknIQxuCn7EAyx0XMK4nhEAgM1HzyGrsEziNERE1BSyCsuw5WgBAGBs5fs+sdBxCa0DvHBj+wAIcemSQyIici5Vk5BvbB+AcD9PidPYDxY6LmJ8L2t1v2RPFgwmTkomInImJrMFP++1FjrszamOhY6LGFS510lBiQHrDpyVOg4RETWivw7l42yRHn5e7hgS3VLqOHaFhY6LUCrktksNF+86JXEaIiJqTIt3XZqE7O7Gj/Z/4k/DhYxNjIBcBmw/fh4nzpVIHYeIiBrB6fNl2Hz0HADgvkQOW12OhY4LaeXrgYEdresqcFIyEZFzWJx0CkIA/TsEItLfS+o4doeFjoupmpT8894zqDCaJU5DRETXQ28y4+c91pWQJ/Rib87VsNBxMQM7BiFUq8bFMiP+zMiTOg4REV2HPzPyUFhqQIhWzZWQa8BCx8Uo5DKMrRzD5aRkIiLH9t1O6/v42J4RXAm5Bnb9U9m8eTNGjRqF0NBQyGQyLFu2rNrjQgjMmjULISEh8PDwwODBg3H06FFpwjqQe3uGQyGXYffJCzhytljqOERE1ACH8oqw++SFyj9guYFnTey60CktLUW3bt3wySefXPXxt99+Gx9++CE+++wz7Nq1C15eXhg2bBgqKiqaOaljaalRY0hn6zoL3P+KiMgxVb1/D+ncEi01aonT2C+7LnRuvfVWzJkzB3fccccVjwkhMG/ePPznP//B7bffjtjYWHz77bfIycm5oueHrnRf5aS1X5PPoMxgkjgNERHVR6nehN+SswEAE26IlDiNfbPrQqc2mZmZyMvLw+DBg233abVa9OrVCzt27KjxeXq9HkVFRdVurqhfuwBE+HmiuMKEVWm5UschIqJ6WJ6agxK9CVEBXujT1l/qOHbNYQudvDzrFUMtW1Zf6rply5a2x65m7ty50Gq1tlt4uGuOa8rlMluvzmKuqUNE5DCEELZJyPclRkAul0mcyL45bKHTUC+88AJ0Op3tlpWVJXUkydzTIwxKhQxpWReRfkYndRwiIqqD1KyLOJBbBHc3Oe7uESZ1HLvnsIVOcHAwAODs2eobVJ49e9b22NWoVCpoNJpqN1fl763C8JgQAMC3O05KG4aIiOrku53WXviRsSFo4eUucRr757CFTlRUFIKDg7FhwwbbfUVFRdi1axd69+4tYTLHMql3awDAirQcXCg1SBuGiIhqdbHMgFX7cgAA43txEnJd2HWhU1JSgtTUVKSmpgKwTkBOTU3F6dOnIZPJMH36dMyZMwcrVqxAeno6Jk2ahNDQUIwePVrS3I6ke4QvuoRqoDdZsGSP6w7jERE5gl/2noHeZEHnEA26R/hKHcch2HWhs2fPHsTHxyM+Ph4AMGPGDMTHx2PWrFkAgGeffRZPPPEEHn74YfTs2RMlJSX4888/oVZzPYG6kslkmFzZq7No5ymYLULaQEREdFUWi8DiyrVzJtwQAZmMk5DrQiaEcOlPtqKiImi1Wuh0Opedr1NuMOOGuRugKzfi68kJGNS55bWfREREzWrTkXOY/E0SfFRu2PHiIHir3KSOJKm6fn7bdY8ONQ8PdwXu7Wm9zH7hDu5/RURkjxZuPwkAuKtHmMsXOfXBQocAABN6RUImAzYfOYfMglKp4xAR0T+cOl+Kvw/nAwAm9eYk5PpgoUMAgAh/T9zUMQgAsIi9OkREduW7nacgBNC/QyDaBHpLHcehsNAhm4mVfyX8vDeL+18REdmJMoMJP+22XhV7fx/25tQXCx2yGdA+EJH+1v2vlqXkSB2HiIhg3deqqMKECD9PDOwQJHUch8NCh2zkchkmVu6C++2Ok3DxC/KIiCQnhLBNQp7UO5L7WjUACx2q5p4e4VAr5TiUV4zdJy9IHYeIyKUlZRbiUF4xPJQK3NPDNTehvl4sdKgaracSo+NaAQAWcv8rIiJJVb0Pj45vBa2nUtowDoqFDl2halLymow8nC2qkDgNEZFrytWVY81+68bVkzkJucFY6NAVuoRqkRDZAiaLwPeVy40TEVHzWrzzNMwWgV5RfugU7Jor9zcGFjp0VZP6tAYAfJ90GnqTWdowREQuRm8y44ck6x+a91e+H1PDsNChq7q1azBaalQ4V6zH6n25UschInIpq/fl4nypASFaNYZEc//B68FCh65KqZBjUuWu5t9sy+Sl5kREzahq38EJN0TCTcGP6uvBnx7VaFxiBFRucmRkF2HPKV5qTkTUHFKzLiIt6yLcFXKM7clLyq8XCx2qkZ+XO+6It15q/s3WTInTEBG5hqr325HdQuDvrZI4jeNjoUO1eqBvFABgzf48ZBWWSZyGiMi55erK8Xu6dV7klMr3X7o+LHSoVh2DfdC3nT8sAli0k7uaExE1pW93nIKp8pLyrq20UsdxCix06Jqq/qr4Iek0SvXc1ZyIqCmUGUy2tcum9mNvTmNhoUPXdFPHILSu3NX8t+QzUschInJKvyVnQ1duRKS/JwZ15iXljYWFDl2TXC6zLVg1f9tJWCy81JyIqDFZLALfbLNOQr6/T2souEt5o2GhQ3Vyd0I4fFRuOFFQik1Hzkkdh4jIqWw6cg4nzpXCR+WGexJ4SXljYqFDdeKtcsOYyvUcqv7qICKixlH1vnpvz3B4q9wkTuNcWOhQnd3fpzXkMmDL0QIcPVssdRwiIqdwOK8YW44WQC4DJnNfq0bHQofqLNzPE4MrJ8jN335S2jBERE6iaoHAW7oGI9zPU+I0zoeFDtXLlMpLHn9LPoMLpQaJ0xARObbzJXosTc0GwAUCmwoLHaqXXlF+iA7RoMJowfdJp6WOQ0Tk0BbvOg2DyYJuYVr0iGwhdRynxEKH6kUmk+Gh/ta/OuZvOwm9ySxxIiIix6Q3mfFt5S7lU/pFQSbjJeVNgYUO1dvI2FAEa9QoKNFjeUqO1HGIiBzSqrRcFJTo0VKjwvCYEKnjOC0WOlRvSoUcU/q1BgB8seUEFxAkIqonIQS+qpyEPKl3aygV/DhuKvzJUoOMTYyAt8oNx/JLuIAgEVE9bT1WgIO5RfBQKjC+V4TUcZwaCx1qEI1aiXGJ1gUEv9xyQuI0RESO5YvN1vfNe3uGw9fTXeI0zo2FDjXYA32j4CaXYfvx88jI1kkdh4jIIezP0dkWCOQu5U2PhQ41WKivB0bEWifQsVeHiKhuvqzszRkeE8IFApsBCx26Lg/d2AYAsGpfLrIvlkuchojIvmVfLMfKfbkAgEf6t5U4jWtgoUPXpWsrLfq09YfZIjB/Kzf7JCKqzTdbM2G2CPRu44+YMK3UcVwCCx26bg/1t/bq/JB0Grpyo8RpiIjsk67ciB8rV5R/eEAbidO4DhY6dN0GdghE+yBvlBrMthcxERFV9/2u0yg1mNGxpQ8GdgiUOo7LYKFD1826LYT1r5P5207CYLJInIiIyL7oTWbM32Yd3n+ofxtu99CMWOhQo7g9LhSBPirkFVVgdTq3hSAi+qflqTnIL9YjWKPGbd1CpY7jUljoUKNQuSlwf5/WAIDPN52AENwWgogIACwWYbuk/IG+reHuxo/e5sSfNjWaCb0i4eWuwKG8Yvx9OF/qOEREdmHjkXwczS+Bt8oN47jdQ7NjoUONRuupxIQbIgEAn248LnEaIiL78Pkma2/OuMRwaNRKidO4HhY61Kim9IuCu0KO3ScvYPfJQqnjEBFJKuX0BezKLISbXIYH+nK7BynUu9ApLy9Hdnb2Fffv37+/UQKRY2upUeOuHmEA2KtDRPS/yvfB0fGtEOrrIXEa11SvQueXX35B+/btMWLECMTGxmLXrl22xyZOnNjo4cgxPdK/DeQy4K9D+TiYWyR1HCIiSRzOK8a6A2chkwH/GsDtHqRSr0Jnzpw52Lt3L1JTUzF//nxMnToV33//PQBIcpWN2WzGyy+/jKioKHh4eKBt27Z47bXXeMWPxFoHeGF4jHWzz882sVeHiFxT1fvfLV2C0S7IW+I0rsutPo2NRiNatmwJAOjRowc2b96MO+64A8eOHZNk8aO33noLn376KRYuXIguXbpgz549eOCBB6DVavHkk082ex665NGBbbFqXy5WpuXg30M6IsKfO/QSkevIKizDijTrmmKPDWwncRrXVq8enaCgIOzbt8/2bz8/P6xbtw4HDx6sdn9z2b59O26//XaMGDECrVu3xt13342hQ4ciKSmp2bNQdV1CtRjQIRAWAXy+mb06RORaPt98HGaLwI3tA7h5p8TqVegsWrQIQUFB1e5zd3fHDz/8gE2bNjVqsLro06cPNmzYgCNHjgAA0tLSsHXrVtx66601Pkev16OoqKjajZrGYwOtY9I/7z2D/OIKidMQETWP/KIKLNlzBgAw7Sb25kitXoVOWFgYgoODq923fv16AEDfvn0bL1UdPf/88xg7diw6deoEpVKJ+Ph4TJ8+HePHj6/xOXPnzoVWq7XdwsPDmzGxa0mM8kP3CF8YTBZ8s/Wk1HGIiJrF11szYTBZ0D3CF72i/KSO4/Kuex2dESNGYMaMGTAYDI2Rp16WLFmCxYsX4/vvv0dycjIWLlyId955BwsXLqzxOS+88AJ0Op3tlpWV1YyJXYtMJrONTS/eeQpFFUaJExERNS1dmRHf7TwFwNqbw807pXfdhc7mzZuxatUqJCQkICMj46ptcnNzcdddd13vt7rCzJkzbb06MTExmDhxIp5++mnMnTu3xueoVCpoNJpqN2o6N3cKQoeW3ijWm7Boxymp4xARNamFO06i1GBGp2Af3Nwp6NpPoCZ33YVOr169kJycjISEBPTs2RPvvfee7TGLxYIDBw5g1qxZ2LJly/V+qyuUlZVBLq9+CAqFAhaLpdG/FzWMXC7Do5VzdeZvy0SF0SxxIiKiplFmMGH+tkwA1itP2ZtjH+p1eXlNvL298e6778LT0xMzZ87EDz/8YCty9Ho9IiMja+1laahRo0bh9ddfR0REBLp06YKUlBS89957mDJlSqN/L2q4kbGheHftEZy5UI4fkk5zGXQicko/JGXhQpkRkf6eGFG5lhhJ77p7dL766itEREQgICAACxYsQGJiItzc3JCSkoIHH3wQhYWFyMzMxNSpUxsjbzUfffQR7r77bjz22GPo3LkznnnmGTzyyCN47bXXGv17UcMpFXLbXJ3PNh1nrw4ROR2DyYIvN1s373ykf1u4KbiVpL2QietcRjgoKAh33XUXpk+fjg4dOti66t5//328+OKLGDduHD7++GN4etrngnFFRUXQarXQ6XScr9OE9CYzBv7fRuTqKvDa6K6YWLnLORGRM/gh6TRe+C0dQT4qbHnuJqjcFFJHcnp1/fy+7pJz4MCBePXVV9GxY8dq45FPP/00kpKSsGfPniv2xSLXo3JT2ObqfPr3MRhMnEdFRM7BaLbgk7+PAQAe7t+GRY6due5CZ8mSJbZtIS4XExOD3bt3Y+TIkejfv//1fitycGMSwhHko0KOrgK/JZ+ROg4RUaNYmpKNMxfKEeDtjvG92Fttb5p8EFGlUmHevHlYtWpVU38rsnNqpQKPVO7g+8nGYzCa2atDRI7NdFlvjoc7e3PsTbPNlhoyZEhzfSuyY/clRiDA2x1ZheVYnpojdRwiouuyPDUHp86Xwc/LHRM499AucVo4NSsPdwUeurENAOCTv4/BxF4dInJQJrMFH1f25jx0Yxt4ujfKii3UyFjoULObcEMkWngqkVlQilX7cqWOQ0TUIKv25SKzoBQtPJWY1Ju9OfaKhQ41Oy+VGx6s7NX5+O9jMFuua4UDIqJmZ7YIfPjXUQDAgze2gZeKvTn2ioUOSWJS70ho1G44ll+CPzLYq0NEjmV1ei5OnCuF1oO9OfaOhQ5JwketxJR+1q0gPv7rGCzs1SEiB2GxCHy0wdqbM7VfFHzUSokTUW1Y6JBkHugTBW+VGw7lFWPN/jyp4xAR1ckfGXk4ml8CH7Ub7u/bWuo4dA0sdEgyWk8lHqh8k5i3/ih7dYjI7lksAh9Vzs2Z0jcKGvbm2D0WOiSpB/u1gY/aDYfPFmN1OufqEJF9W3sgD4fyiuGjcsOUvlFSx6E6YKFDktJ6KvFgP+sVWPPWH+EVWERkt8wWgffWHQEAPNC3NbSe7M1xBCx0SHJT+rWGr6cSx8+VYnlqttRxiIiuatW+HBw5WwKN2g1TK5fIIPvHQock56NW4uH+1jeNDzYc5WrJRGR3TGYLPlhvnZvzcP820HqwN8dRsNAhuzC5d2v4e7nj1Pky/JbMXh0isi9LU7JxonIV5Ps5N8ehsNAhu+ClcsO/Knc2/2DDURhM7NUhIvtgMFnwQeW6OY8ObAtvroLsUFjokN2YcEMkAn1UyL5Yjp/3Zkkdh4gIAPDz3iycuVCOQB8VJt7QWuo4VE8sdMhueLgr8NhAa6/Ox38dQ4XRLHEiInJ1FUYzPtpg3aF82sC28HBXSJyI6ouFDtmVcYkRCNGqkaurwI9Jp6WOQ0Qu7oek08grqkCIVo2xiRFSx6EGYKFDdkWtVGDaTe0AAJ9sPM5eHSKSTLnBjE/+Pg4AeOLm9lAr2ZvjiFjokN0ZkxCOVr4eOFesx7c7Tkodh4hc1Lc7TqKgRI9wPw/ckxAmdRxqIBY6ZHfc3eR4anB7AMD/Nh5HUYVR4kRE5GpK9CZ8tsnam/PUoA5QKvhx6ah45sgu3RnfCu2CvHGxzIgvNp2QOg4RuZivtpzAhTIj2gR6YXRcqNRx6Dqw0CG75KaQ45mhHQEAX2/NRH5xhcSJiMhVFJTo8eVm6x9YM4Z0gBt7cxwazx7ZrWFdWiIu3BflRjM+/uuY1HGIyEV8/NcxlBrMiGmlxfCuIVLHoevEQofslkwmw3O3dAIAfL/rNE6fL5M4ERE5u6zCMizedQoA8NwtnSCXyyRORNeLhQ7Ztd5t/dG/QyBMFoF31x2WOg4RObl31x6G0SxwY/sA9GsfIHUcagQsdMjuPTvMOldneWoO9ufoJE5DRM7qQE4RlqflAICtN5kcHwsdsntdW2kxqpv1qod31rBXh4iaxttrDkEIYGRsCLq20kodhxoJCx1yCP8e0gFuchn+PnwOu06clzoOETmZHcfPY+Phc3CTy2xXfJJzYKFDDqF1gBfu7RkOAHh7zWEIISRORETOQgiBN/88BMC6317rAC+JE1FjYqFDDuOpQe2hVsqx99QFrDtwVuo4ROQk1uzPQ1rWRXgoFXhiUDup41AjY6FDDiNIo8aUvlEAgDf/PASj2SJxIiJydCazBW9Xzv178MYoBPmoJU5EjY2FDjmURwe2hb+XO06cK8WPSaeljkNEDm7JnjM4ca4ULTyVeLh/G6njUBNgoUMOxUetxPTKDT/fX3+UG34SUYOV6E14r3J9ridubg8ftVLiRNQUWOiQwxmbGIE2gV4oLDXgs43HpY5DRA7q043HUFBiQFSAFybcECl1HGoiLHTI4SgVcrx4a2cA1g0/sy+WS5yIiBxNzsVyfLUlEwDw/K2d4O7Gj0NnxTNLDmlQ5yDc0MYPepOFiwgSUb3935rD0JssSIzyw9DollLHoSbEQocckkwmw0vDowEAS1OykX6GW0MQUd3sO3MRS1OyAQD/GdEZMhk37nRmLHTIYcWEaXFnfCsAwJzVB7iIIBFdkxACc1YfBADcEd8KsWG+0gaiJsdChxzav4d1hMpNjl2ZhVh/MF/qOERk59YeOIukzEKo3OSYOYxbPbgCFjrk0Fr5emBqP+signP/OMhFBImoRgaTBW/+Yd3q4cEboxDq6yFxImoOLHTI4f1zEcHFO09JHYeI7NR3O08hs6AUAd7ueHQgt3pwFSx0yOH5qJV4ekgHANZFBC+UGiRORET25mKZAR/+dRQAMGNIR3ir3CRORM3F4Qud7OxsTJgwAf7+/vDw8EBMTAz27NkjdSxqZuMSI9Ap2Ae6ciPeXcfLzYmouvfXHcHFMiM6tPTGmIQwqeNQM3LoQufChQvo27cvlEol/vjjDxw4cADvvvsuWrRoIXU0amYKuQyv3tYFAPD9rtM4mFskcSIisheH8orw3S7r3nivjOoCN4VDf/RRPTl0391bb72F8PBwzJ8/33ZfVFSUhIlISje08ceImBCsTs/F7JX78cNDN3B9DCIXJ4TA7BUHYLYI3NIlGH3bBUgdiZqZQ5e1K1asQEJCAu655x4EBQUhPj4eX375Za3P0ev1KCoqqnYj5/HC8E5Qucmx80Qh/sjIkzoOEUnsz4w87DhxHio3OV4a0VnqOCQBhy50Tpw4gU8//RTt27fHmjVr8Oijj+LJJ5/EwoULa3zO3LlzodVqbbfw8PBmTExNLayFJx4Z0BYA8Prqg6gwmiVORERSqTCabYsDPtK/DcL9PCVORFKQCQdeTtbd3R0JCQnYvn277b4nn3wSu3fvxo4dO676HL1eD71eb/t3UVERwsPDodPpoNFomjwzNb1ygxmD3t2IHF0FZgzpgCcHtZc6EhFJ4IP1R/H++iMI1aqx4d8D4eGukDoSNaKioiJotdprfn47dI9OSEgIoqOjq93XuXNnnD59usbnqFQqaDSaajdyLh7uCrww3NpF/b+Nx5DD3c2JXE72xXJ8uukYAOCF4Z1Z5Lgwhy50+vbti8OHq19KfOTIEURGRkqUiOzFyNgQJLb2Q4XRgrmVK6ESket44/eDqDBadycfGRsidRySkEMXOk8//TR27tyJN954A8eOHcP333+PL774AtOmTZM6GklMJpNh1qhoyGTAyrQc7DxxXupIRNRMdhw/j9X7ciGXAa+O6sKrL12cQxc6PXv2xNKlS/HDDz+ga9eueO211zBv3jyMHz9e6mhkB7q20mJcYgQAYNbyDO6DReQCTGYLZq/cDwC4r1cEokM5PcHVOfQ6OgAwcuRIjBw5UuoYZKeeHdYRf2bk4cjZEszflomH+7eVOhIRNaEF20/iUF4xfD2VmDGEu5OTg/foEF2Lr6c7nr+1EwBg3vqjyNVxYjKRs8rVleP9dUcAAM/f0gl+Xu4SJyJ7wEKHnN7d3cOQENkCZQYzXlt1QOo4RNRE5qw6iFKDGd0jfDEmgWukkRULHXJ6crkMr43uCoVcht/T87DxcL7UkYiokW06cg6r060TkOeMjoFczgnIZMVCh1xC5xAN7u/TGgDwyor9XDGZyIlUGM14ZXkGAOD+PlGcgEzVsNAhlzF9cHu01Khw6nwZPt90Quo4RNRIPtt0HCfPl6GlRoWnh3AldKqOhQ65DB+1Ev8ZYV1J+5ONx3DqfKnEiYjoep0sKMX/Nh4HALw8Mho+aqXEicjesNAhlzIyNgT92gXAYLLglRX74cBbvRG5PCEEZq3YD4PJghvbB2BEDFdApiux0CGXIpPJMPv2LnBXyLHx8Dms3JcrdSQiaqDf0/Ow+cg5uCvkmH0bV0Cmq2OhQy6nbaA3HrvJunDgf1fux8Uyg8SJiKi+dGVGvLLCugLyvwa0QZtAb4kTkb1ioUMu6dGBbdEuyBsFJQa8vvqg1HGIqJ7m/nEQBSV6tAn0wmM3tZM6DtkxFjrkklRuCrx5ZwwA4Oe9Z7D9WIHEiYiorrYfL8CPu7MAAG/eGQu1UiFxIrJnLHTIZSW09sOEG6ybfr6wNJ1r6xA5gAqjGS/+lg4AGN8rAolRfhInInvHQodc2rO3dLKtrfPBhqNSxyGia/hgw1HbmjnPVe5jR1QbFjrk0jRqJf57e1cAwBebT+BATpHEiYioJvtzdPhis3Wxz//e3hUarplDdcBCh1zesC7BuKVLMMwWgRd+2wezhWvrENkbk9mC539Nh9kicGvXYAzrEix1JHIQLHSIAMy+vQt81G5IO6PD/G2ZUschosss2H4S6dk6+KjdMPu2LlLHIQfCQocIQEuNGi/c2hkA8M7aw8gs4PYQRPbiZEEp3ll7GADw0vDOCNKoJU5EjoSFDlGlcYnh6NvOHxVGC2b+nMYhLCI7YLEIzPwlDRVGC3q38ce9PcOljkQOhoUOUSWZTIa37oqFl7sCe05d4BAWkR2Yv/0kdp+8AE93Bd6+O5bbPFC9sdAh+oewFp54qXKH8/9bcxgnzpVInIjIdZ04V4K3/zwEAHhxeGeE+3lKnIgcEQsdosuMSwzHje0DoDdZMPMXXoVFJAWzRWDmL/ugN1nQr10AxveKkDoSOSgWOkSXkclkePOuWHir3LCXQ1hEkvhmayb2nroAb5Ub3rwrhkNW1GAsdIiuopWvB14aYb0K6//WHMZxDmERNZtj+SX4v8qrrP4zojPCWnDIihqOhQ5RDcb2/McQFq/CImoWZovAMz+nwWCyoH+HQF5lRdeNhQ5RDf45hJV8+iI+23Rc6khETu/zzceRmnURPio3vHknh6zo+rHQIapFK18PvFq5Cuv7644g/YxO4kREzisjW4f31h4BALw8Khqhvh4SJyJnwEKH6Bru6t4Kw2OCYbIIPPVTCsoNZqkjETmdcoMZT/2YApNF4JYuwbinR5jUkchJsNAhugaZTIbXR8egpUaFE+dK8cbvB6WOROR03vj9II6fK0WQjwpzOWRFjYiFDlEdtPByxzv3dAMALNp5Cn8fypc4EZHz+PtQPhbtPAUAeOeebmjh5S5xInImLHSI6ujG9oGY0jcKADDzlzQUlOglTkTk+ApK9Jj5SxoA4IG+rdG/Q6DEicjZsNAhqodnb+mIDi29UVBiwPO/pkMIXnJO1FBCCDz/6z4UlBjQoaU3nrulk9SRyAmx0CGqB7VSgXn3xsNdIcf6g2fxfdJpqSMROawfkrKw/mA+3BVyzLs3HmqlQupI5IRY6BDVU3SoBs8M6wAA+O/KAzicVyxxIiLHc+RsMf67aj8AYOawjogO1UiciJwVCx2iBniwXxv07xAIvcmCx79PRpnBJHUkIodRbjBj2uJkVBgtuLF9AKb2i5I6EjkxFjpEDSCXy/DemG4I8lHhaH4JXl2xX+pIRA5j9sr9OJpfgkAfFd6/Nw5yOS8lp6bDQoeogQK8VZg3Ng4yGbBkzxksT82WOhKR3Vuemo0fd2dBJgM+uDcOAd4qqSORk2OhQ3Qd+rQNwBM3twcAvPhbOjILSiVORGS/MgtK8eJv6QCAJ25ujz7tAiRORK6AhQ7RdXpqUHv0ivJDqcGMx79Pht7ELSKILqc3WV8fpQYzekX54alB7aWORC6ChQ7RdVLIZfhgbDz8vNyxP6cIb6zmFhFEl5v7+yHszymCn5c7PhgbDwXn5VAzYaFD1AiCtWq8W7lFxMIdp7AiLUfiRET2Y9W+HCzYfhIA8O493RCsVUsbiFwKCx2iRnJTpyBMu6ktAOC5X/bhyFmur0N05Gwxnv1lHwDgXwPa4qZOQRInIlfDQoeoEc0Y0hH92gWg3GjGvxbtRVGFUepIRJIprjDiX4v2osxgRt92/nhmaAepI5ELYqFD1Iis83XiEKpV40RBKZ5Zksb9sMglCSHwzM9pOFFQihCtGh+OjYebgh851Pz4W0fUyPy9VfjfhB5wV8ix9sBZfLbphNSRiJrd55tPYM3+s3BXyPHphB7w53o5JBEWOkRNIC7cF6/e1gUA8H9rDmH7sQKJExE1n+3HCvD2n4cAAK/cFo24cF9pA5FLc6pC580334RMJsP06dOljkKEcYnhuLtHGCwCePyHFGRfLJc6ElGTy75Yjid+SIFFAHf3CMN9iRFSRyIX5zSFzu7du/H5558jNjZW6ihEAACZTIY5o7uiS6gGhaUGPLhwDzf/JKdWZjDhoYV7cL7UgOgQDeaM7gqZjOvlkLScotApKSnB+PHj8eWXX6JFixZSxyGyUSsV+GJSAgK83XEwtwj/XpIGi4WTk8n5WCzWyccHcovg7+WOLyb1gFqpkDoWkXMUOtOmTcOIESMwePDga7bV6/UoKiqqdiNqSq18PfDZhB5QKmT4IyMPH2w4KnUkokb34V9H8Xt6HpQKGT6f2ANhLTyljkQEwAkKnR9//BHJycmYO3dundrPnTsXWq3WdgsPD2/ihERAQms/vH5HDADggw1HsXpfrsSJiBrP7+m5mLfeWsC/PjoGCa39JE5EdIlDFzpZWVl46qmnsHjxYqjVdVtS/IUXXoBOp7PdsrKymjglkdWYhHBM7RcFAPj3z6nIyNZJnIjo+mVk6zBjSSoAYErfKIzpyT8eyb7IhAOvZrZs2TLccccdUCgujQObzWbIZDLI5XLo9fpqj11NUVERtFotdDodNBpNU0cmF2cyWzBl4R5sPnIOIVo1lj/eF0E+3PeHHFN+cQVGf7wNOboK9O8QiG8mJ3BRQGo2df38dujfyEGDBiE9PR2pqam2W0JCAsaPH4/U1NRrFjlEzc1NIcdH4+LRJtALuboKTF3AK7HIMZUZTHhw4R7k6CrQJsALH43jysdknxz6t9LHxwddu3atdvPy8oK/vz+6du0qdTyiq9J6KPHN5J5o4alEerYOT3yfApPZInUsojozWwSe/CEF+87o0MJTia/v7wmth1LqWERX5dCFDpGjah3gha8m94TKTY4Nh/Ixe+UB7olFDkEIgVdX7Mf6g/lQucnx1eSeiArwkjoWUY3cpA7Q2DZu3Ch1BKI66RHZAvPujcNj3ydj0c5TCPfzwMP920odi6hWX245gUU7T0EmA+bdG4cekVy7jOwbe3SIJHRrTAheGt4ZAPDG74ewal+OxImIarZ6Xy7e+N26h9VLwzvj1pgQiRMRXRsLHSKJTe0Xhfv7tAYAzFiShqTMQmkDEV1FUmYhnq68jPz+Pq1tSyUQ2TsWOkQSk8lkeHlkNIZGt4TBZMHUBbu5xg7ZlYxsHaYu2A2DyYKh0S3x8sho7mFFDoOFDpEdUMhl+HBcPBKj/FCsN+H++UnILCiVOhYRMgtKcf/8JBTrTUiM8sOH4+KhkLPIIcfBQofITqiVCnw1OQFdQjUoKDFgwle7kKsrlzoWubBcXTkmfLULBSUGdAnV4KvJCdyokxwOCx0iO6JRK7FwSiLaBHgh+2I5Jn6dhMJSg9SxyAUVlhow8eskZF8sR5sALyyckgiNmmvlkONhoUNkZwK8Vfh2aiJCtGocyy/B/fOTUKLn6snUfEr0JjwwPwnH8ksQolXj26mJCPBWSR2LqEFY6BDZobAWnlg0NREtPJXYd0aHKfN3o5TFDjWD0soiJ61y1eNFUxMR1sJT6lhEDcZCh8hOtQvywbdTesFH7Yakk4WYsmA398WiJlVmMOGBBbux++QF+KjdsHBKItoF+Ugdi+i6sNAhsmMxYVosmtoLPio37MosxNQFe1BuMEsdi5xQucGMqQv2ICmzED4qNyya2guxYb5SxyK6bix0iOxcXLgvFk5NhLfKDTtOnMdD3+5BhZHFDjWeCqMZD327BztOnIe3yg0LpyYiLtxX6lhEjYKFDpED6B7RAgun9ISXuwJbjxWw2KFGU2E045FFe7H1WAE83RVY8EBPdI/g/lXkPFjoEDmIHpF+WDAlEZ7uCmw5WoApCzhBma6PdeLxbmw6cg4eSgUWPJCIhNZ+UscialQsdIgcSM/Wfph/v7VnZ/vx85jw9S7oyoxSxyIHpCszYsLXu2zDVQse6InEKBY55HxY6BA5mF5t/LH4oRug9VAi5fRFjP1yJwpK9FLHIgdyvkSPcV/uRMrpi9B6KLH4wV7o1cZf6lhETYKFDpEDigv3xU+P3IAAbxUO5hZhzOc7uF0E1UmergJjPt+BA7lFCPB2x48P34BunHhMToyFDpGD6hSswZJHbkCoVo0T50pxz2c7cOJcidSxyI4dP1eCuz/bjuPnShGiVWPJI73ROUQjdSyiJsVCh8iBtQn0xpJ/9UZrf0+cuVCOuz7djuTTF6SORXZo76kLuPvT7ThzoRyR/p5Y8khvtAn0ljoWUZNjoUPk4MJaeOLnf/VBbJgWF8qMGPfFTqzdnyd1LLIja/fn4b4vd+JCmRHdwrT49dE+CPfjtg7kGmRCCCF1CCkVFRVBq9VCl5MDjeYqXbgKBaBWX/p3aWnNX0wuBzw8Gta2rAyo6VTIZICnZ8PalpcDFkvNOby8Gta2ogIw17KOS33aenpacwOAXg+Yarlkuj5tPTysP2cAMBgAYy1XJ9WnrVpt/b2ob1uj0dq+JioV4OZW/7YmE6DXo8xgwowladh0+BzkMuA/IzpjXK9IwN0dUCqrta3RP9uazdZzVxOl0tq+vm0tFuvvWmO0dXOz/iwA62uirKxx2tbndW/n7xE/JJ3GnFUHYBFA/w4BeH9sPDx9//Fex/cI6/+7wHtEjRz0PcL2+a3TXf3zu4pwcTqdTgAQOuvbwpW34cOrP8HT8+rtACEGDKjeNiCg5rYJCdXbRkbW3DY6unrb6Oia20ZGVm+bkFBz24CA6m0HDKi5radn9bbDh9fc9vJfq7vvrr1tScmltpMn1942P/9S28ceq71tZualts88U3vbjIxLbV95pfa2SUmX2r79du1t//77UtuPP6697apVl9rOn1972yVLLrVdsqTWtuavv7nUdtWq2r/uxx9favv337W3ffvtS22Tkmpv+8orl9pmZNTe9plnLrXNzKy97WOPXWqbn19728mTL7UtKam97d13i2pqa8v3COuN7xGXbg72HiHmz7/U1oHeI2yf3zqdqA2Hroic3MIdJ1HChQVdillInYDIfnDoikNXDWvLbun6t23mbukVadl4edl+lAo52oT64qtJPRGhdXeYbuk6teXQldU/XvdZhWWY9uVWHD1bDHc3OV6/oytGxoZetS0Avke48HuEjZMPXbHQqesYH5EDSs26iIe/3YP8Yj18PZX43/ju6NM2QOpY1ER2njiPR7/biwtlRgT5qPDFpARuzklOq66f3xy6InJiceG+WPF4P8SGaXGxzIiJXydhwbZMuPjfN05HCIEF2zIx4atduFBmRGyYFise78cihwgsdIicXnDlwnC3x4XCbBF4deUBPP5DCooruEeWMyiuMOLx71Pw6soDMFkERnULxU8P90awVn3tJxO5ABY6RC5ArVRg3r1xeHlkNNzkMqzel4vbPt6Gg7lFUkej63Aorwi3f7wNq9Nz4SaX4ZVR0fhwbBw83BVSRyOyGyx0iFyETCbD1H5R+OmR3gjRqpFZUIrRn2zDkt1ZUkejehJC4Je9ZzD6k204UWDdzuGnR3rjgb5RkFVNxCUiACx0iFxOj8gWWP3kjRjQIRB6kwXP/roPT/yQAl0Zh7IcwcUyAx7/IQXP/JyGCqMFN7YPwOonb0SPyBZSRyOyS7zqilddkYuyWAT+t/EY3l9/FGaLQLBGjXfHdEPfdrwqy15tP1aAGUvSkFdUATe5DNMHt8ejA9tBIWcvDrkeXl5eRyx0yNWlnL6AGUvSkFlgXdNlSt8oPHtLR6iVnOdhLyqMZry37gi+2HwCABAV4IV598ahG6+qIhfGQqeOWOgQAWUGE15ffRCLd50GALQL8sZbd8WgR6SfxMlo98lCPPfrPpw4Zy1ExyVG4OWRneHp7iZxMiJpsdCpIxY6RJf8degsnv0lHQUleshkwIRekXj2lo7wUSuljuZySvQmvP3nIXy74xQAINBHhddHd8XQLsESJyOyDyx06oiFDlF1F8sMeH31Qfy89wwAIFijxmuju2JIdEuJk7mODQfP4uVlGcjRWZfMH5MQhpeGR0PryYKTqAoLnTpioUN0dduOFeDFpek4dd66J9TgzkF4aUQ0ogK8rvFMaqjMglL8d+V+/H34HAAgws8Tc++M4QRxoqtgoVNHLHSIalZuMOODDUfx1ZYTMFkElAoZpvSLwuM3teNwViMq1Zvw0V/H8PXWEzCaK3/OfaPw1OD2nItDVAMWOnXEQofo2o7lF+O/qw5i8xFrT0OgjwrPDO2Au7qHwU3B5bgaymi24Oc9ZzBv/RHkF1t3jB7QIRCzRkWjbaC3xOmI7BsLnTpioUNUN0II/HUoH6+tOoCTlcNZbQK88PSQDhgREwI513KpM4tF4PeMXLy79ojtsv4IP0+8PDIagzsHcXVjojpgoVNHLHSI6kdvMmPRjlP45O9juFC5mnLnEA3+PaQDBvFDulYWi7VYnLfhCDKyrfuM+Xm54/Gb2mH8DRFQuXHtIqK6YqFTRyx0iBqmuMKIb7aexFdbTqBYbwIAdAr2wSMD2mBkbCiUHNKyMZktWLUvF59uPI7DZ4sBAF7uCjzcvy2m3hgFbxXn4RDVFwudOmKhQ3R9LpQa8PnmE1i04yRKDWYAQCtfD0ztF4UxPcNd+kO8uMKI35Kz8dXWE8gqLAcAeKvcMOGGSDx0YxT8vVUSJyRyXCx06oiFDlHj0JUZ8d2uU5i/LRMFJQYA1l6L0fGtMOGGSHQOcZ3X1+G8Yny74ySWpmSjrLL48/Nyx9R+UZhwQyS0Hrxijeh6sdCpIxY6RI2rwmjGr8ln8PXWTNu2BQDQPcIXY3tGYFjXYKf8oNeVGbE6PRe/JZ/BnlMXbPe3C/LGpN6RuKdHODzcOQeHqLGw0KkjFjpETUMIgR0nzmPxztNYsz8PJov1rcZdIcdNnQJxe1wr3NwpyKE3Dy03mLHpyDksS8nGX4fyYTBbAAAKuQxDo1tiYu9I9G7jzwnaRE2AhU4dsdAhanr5xRX4ec8ZLEvJxtH8Etv9HkoF+rYLwKDOQbipYxCCtWoJU9bNuWI9/jp0FusOnMWWowXQmyy2xzoF++CO+Fa4Pa6VQxwLkSNziUJn7ty5+O2333Do0CF4eHigT58+eOutt9CxY8c6fw0WOkTNRwiBQ3nFWJGWgxWpOci+WF7t8c4hGtzQxg+Jrf3QM8oPAXYwWfdCqQG7Mgux88R57DxxHofyiqs93srXA8NjgnFHfBiiQ/keQtRcXKLQueWWWzB27Fj07NkTJpMJL774IjIyMnDgwAF4edVtPx4WOkTSEEJgf04R/j6Ujw2H8pF25iIufzdqE+iFrqFadA7RoHOID6JDNAj0UTXJUJAQAvnFehw9W4KMHB0ysnXYn1NkW9Dvn2LDtBjSuSUGR7dEp2AfDk0RScAlCp3LnTt3DkFBQdi0aRP69+9fp+ew0CGyDwUlemw/fh67MwuRlFloW2/mcl7uCoT7eSKshSci/DwR6KNCC08lfD3d0cJTCS+VG9wUMrjJ5XBXyGERAhUmM/RGCyqMZhRXmFBQosf5UgPOFeuRp6vAyfOlOHW+DOVG81W/Z/sgb9zQxh83tPFHrzb20dNE5Orq+vntVAtc6HQ6AICfn1+NbfR6PfR6ve3fRUVFTZ6LiK4twFuF27qF4rZuoQCAi2UGpJy+iAO5RThYecssKEWpwYxDecVXDCE1BrkMCPfzRNdQLbq00qBrqBZdW2nh5+Xe6N+LiJqH0/ToWCwW3Hbbbbh48SK2bt1aY7tXX30Vs2fPvuJ+9ugQ2b8KoxnZF8uRVViGrAvlOFNYhoISAy6WGVBYZsDFMiPKDWaYLBYYTBYYzQIKuQxqpRwqNwVUSjm8VW4I8FbB38sd/t4qtNSo0NrfC5H+1l4idzeu6EzkCFxu6OrRRx/FH3/8ga1btyIsLKzGdlfr0QkPD2ehQ0RE5EBcaujq8ccfx6pVq7B58+ZaixwAUKlUUKk4vk5EROQKHLrQEULgiSeewNKlS7Fx40ZERUVJHYmIiIjsiEMXOtOmTcP333+P5cuXw8fHB3l5eQAArVYLDw8PidMRERGR1Bx6jk5Na1fMnz8f999/f52+Bi8vJyIicjwuMUfHgWs0IiIiaga8jpKIiIicFgsdIiIiclosdIiIiMhpsdAhIiIip8VCh4iIiJwWCx0iIiJyWix0iIiIyGmx0CEiIiKnxUKHiIiInJZDr4zcGKpWVy4qKpI4CREREdVV1ef2tXZJcPlCp7i4GAAQHh4ucRIiIiKqr+LiYmi12hofd+hNPRuDxWJBTk4OfHx8atwktCGKiooQHh6OrKwsp9ws1NmPD3D+Y3T24wOc/xh5fI7P2Y+xKY9PCIHi4mKEhoZCLq95Jo7L9+jI5XKEhYU12dfXaDRO+ctbxdmPD3D+Y3T24wOc/xh5fI7P2Y+xqY6vtp6cKpyMTERERE6LhQ4RERE5LRY6TUSlUuGVV16BSqWSOkqTcPbjA5z/GJ39+ADnP0Yen+Nz9mO0h+Nz+cnIRERE5LzYo0NEREROi4UOEREROS0WOkREROS0WOgQERGR02Khcx0++eQTtG7dGmq1Gr169UJSUlKt7X/++Wd06tQJarUaMTEx+P3335spacPU5/gWLFgAmUxW7aZWq5sxbf1s3rwZo0aNQmhoKGQyGZYtW3bN52zcuBHdu3eHSqVCu3btsGDBgibPeT3qe4wbN2684hzKZDLk5eU1T+B6mjt3Lnr27AkfHx8EBQVh9OjROHz48DWf5yivw4YcnyO9Dj/99FPExsbaFpLr3bs3/vjjj1qf4yjnrkp9j9GRzt/VvPnmm5DJZJg+fXqt7Zr7PLLQaaCffvoJM2bMwCuvvILk5GR069YNw4YNQ35+/lXbb9++HePGjcPUqVORkpKC0aNHY/To0cjIyGjm5HVT3+MDrCtf5ubm2m6nTp1qxsT1U1paim7duuGTTz6pU/vMzEyMGDECN910E1JTUzF9+nQ8+OCDWLNmTRMnbbj6HmOVw4cPVzuPQUFBTZTw+mzatAnTpk3Dzp07sW7dOhiNRgwdOhSlpaU1PseRXocNOT7AcV6HYWFhePPNN7F3717s2bMHN998M26//Xbs37//qu0d6dxVqe8xAo5z/i63e/dufP7554iNja21nSTnUVCDJCYmimnTptn+bTabRWhoqJg7d+5V248ZM0aMGDGi2n29evUSjzzySJPmbKj6Ht/8+fOFVqttpnSNC4BYunRprW2effZZ0aVLl2r33XvvvWLYsGFNmKzx1OUY//77bwFAXLhwoVkyNbb8/HwBQGzatKnGNo72OvynuhyfI78OhRCiRYsW4quvvrrqY4587v6ptmN01PNXXFws2rdvL9atWycGDBggnnrqqRrbSnEe2aPTAAaDAXv37sXgwYNt98nlcgwePBg7duy46nN27NhRrT0ADBs2rMb2UmrI8QFASUkJIiMjER4efs2/WhyNI52/6xUXF4eQkBAMGTIE27ZtkzpOnel0OgCAn59fjW0c+TzW5fgAx3wdms1m/PjjjygtLUXv3r2v2saRzx1Qt2MEHPP8TZs2DSNGjLji/FyNFOeRhU4DFBQUwGw2o2XLltXub9myZY3zGfLy8urVXkoNOb6OHTvim2++wfLly/Hdd9/BYrGgT58+OHPmTHNEbnI1nb+ioiKUl5dLlKpxhYSE4LPPPsOvv/6KX3/9FeHh4Rg4cCCSk5OljnZNFosF06dPR9++fdG1a9ca2znS6/Cf6np8jvY6TE9Ph7e3N1QqFf71r39h6dKliI6OvmpbRz139TlGRzt/APDjjz8iOTkZc+fOrVN7Kc6jy+9eTo2jd+/e1f5K6dOnDzp37ozPP/8cr732moTJqK46duyIjh072v7dp08fHD9+HO+//z4WLVokYbJrmzZtGjIyMrB161apozSJuh6fo70OO3bsiNTUVOh0Ovzyyy+YPHkyNm3aVGMh4Ijqc4yOdv6ysrLw1FNPYd26dXY9aZqFTgMEBARAoVDg7Nmz1e4/e/YsgoODr/qc4ODgerWXUkOO73JKpRLx8fE4duxYU0RsdjWdP41GAw8PD4lSNb3ExES7Lx4ef/xxrFq1Cps3b0ZYWFitbR3pdVilPsd3OXt/Hbq7u6Ndu3YAgB49emD37t344IMP8Pnnn1/R1hHPHVC/Y7ycvZ+/vXv3Ij8/H927d7fdZzabsXnzZnz88cfQ6/VQKBTVniPFeeTQVQO4u7ujR48e2LBhg+0+i8WCDRs21Dj22rt372rtAWDdunW1jtVKpSHHdzmz2Yz09HSEhIQ0Vcxm5UjnrzGlpqba7TkUQuDxxx/H0qVL8ddffyEqKuqaz3Gk89iQ47uco70OLRYL9Hr9VR9zpHNXm9qO8XL2fv4GDRqE9PR0pKam2m4JCQkYP348UlNTryhyAInOY5NNc3ZyP/74o1CpVGLBggXiwIED4uGHHxa+vr4iLy9PCCHExIkTxfPPP29rv23bNuHm5ibeeecdcfDgQfHKK68IpVIp0tPTpTqEWtX3+GbPni3WrFkjjh8/Lvbu3SvGjh0r1Gq12L9/v1SHUKvi4mKRkpIiUlJSBADx3nvviZSUFHHq1CkhhBDPP/+8mDhxoq39iRMnhKenp5g5c6Y4ePCg+OSTT4RCoRB//vmnVIdwTfU9xvfff18sW7ZMHD16VKSnp4unnnpKyOVysX79eqkOoVaPPvqo0Gq1YuPGjSI3N9d2Kysrs7Vx5NdhQ47PkV6Hzz//vNi0aZPIzMwU+/btE88//7yQyWRi7dq1QgjHPndV6nuMjnT+anL5VVf2cB5Z6FyHjz76SERERAh3d3eRmJgodu7caXtswIABYvLkydXaL1myRHTo0EG4u7uLLl26iNWrVzdz4vqpz/FNnz7d1rZly5Zi+PDhIjk5WYLUdVN1KfXlt6pjmjx5shgwYMAVz4mLixPu7u6iTZs2Yv78+c2euz7qe4xvvfWWaNu2rVCr1cLPz08MHDhQ/PXXX9KEr4OrHRuAaufFkV+HDTk+R3odTpkyRURGRgp3d3cRGBgoBg0aZCsAhHDsc1elvsfoSOevJpcXOvZwHmVCCNF0/UVERERE0uEcHSIiInJaLHSIiIjIabHQISIiIqfFQoeIiIicFgsdIiIiclosdIiIiMhpsdAhIiIip8VCh4iIiJwWCx0iIiJyWix0iIiIyGmx0CEip7J161YolUpUVFTY7jt58iRkMhlOnTolYTIikgILHSJyKqmpqejcuTPUarXtvpSUFLRo0QKRkZESJiMiKbDQISKnkpaWhvj4+Gr3paamolu3bhIlIiIpsdAhIqeSmpqKuLi4avelpKRccR8RuQYWOkTkNMxmMzIyMq7o0UlOTmahQ+SiWOgQkdM4fPgwKioqEBoaartvx44dyM7OZqFD5KJY6BCR00hNTQUAfPTRRzh69Cj++OMPTJo0CQBgMBgkTEZEUmGhQ0ROIzU1FcOGDcOJEycQExODl156CbNnz4ZGo8GHH34odTwikoBMCCGkDkFE1BiGDRuGnj17Ys6cOVJHISI7wR4dInIaaWlpiImJkToGEdkRFjpE5BTy8vJw9uxZFjpEVA2HroiIiMhpsUeHiIiInBYLHSIiInJaLHSIiIjIabHQISIiIqfFQoeIiIicFgsdIiIiclosdIiIiMhpsdAhIiIip8VCh4iIiJwWCx0iIiJyWv8PonpJFq4Mf8cAAAAASUVORK5CYII=" }, "metadata": {}, "output_type": "display_data" }, { "name": "stdout", "output_type": "stream", "text": [ "Calculated mean from LS: 2.000±0.577\n" ] } ], "execution_count": 2 }, { "metadata": {}, "cell_type": "markdown", "source": "The results match those obtained using both the analytical solution and the MCMC algorithm. For this simple problem, the LS method is significantly faster than MCMC. However, MCMC is typically more effective than LS for complex or nonlinear models, non-Gaussian errors, multimodal distributions, or when prior knowledge needs to be included. In contrast, LS is quicker and easier for linear models with Gaussian errors but falls short when dealing with more complicated scenarios.", "id": "536e3a737bd725cf" }, { "metadata": {}, "cell_type": "markdown", "source": [ "## Explanation of affine-invariant algorithms (used in emcee)\n", "### Affine invariance\n", "While the Metropolis-Hasting implementation may be particularly useful for some distributions, it fails for badly scaled or highly skewed distributions. In 2010, Goodman and Weare suggested the use of affine invariant many particle MCMC samplers, where they make use of affine transformations of the form \n", "$$\n", " y \\longmapsto Ax +b. \\tag{6.2}\n", "$$\n", "Furthermore, if a sample $X$ has probability density $p(x)$, then the affine transformed $Y = AX + b$ has probability density $p(y) = p(Ax+b) \\propto p(x)$.\n", "\n", "Rescaling the distribution allows for easier computations with less customization. For example, the Metropolis MCMC application would have difficulty analyzing an anisotropic density like\n", "\n", "\n", "$$\n", " p(x) \\propto \\exp{\\left(-\\frac{(x_1 - x_2)^2}{2a} - \\frac{(x_1 + x_2)^2}{2}\\right)}\n", "$$\n", "\n", "as it will have to take small steps of the order of $\\sqrt{a}$ in all directions due to its anisotropic nature. Ideally, it would like to take bigger steps in the stretched out direction to increase the efficiency, but it is limited by these small steps. This results in a long burn-in time until the Markov chain reaches equilibrium.\n", "\n", "However, using the suited affine transformation \n", "\n", "$$\n", " y_1 = \\frac{x_1 - x_2}{\\sqrt{a}}, \\hspace{1cm} y_2 = x_1 + x_2\n", "$$\n", "\n", "we can alter this problem into a well-scaled case\n", "\n", "\n", "$$\n", " p(y) \\propto \\exp{\\left(\\frac{-(y_1 + y_2)^2}{2}\\right)}. \\tag{6.3}\n", "$$\n", "Where $p(y)$ is much easier to sample than the initial probability density of $p(x)$. \n", "\n", "For an affine invariant algorithm, the two cases of $p(x)$ and $p(y)$ are the same and equally difficult, such that it does not depend on the parameter $a$ that determines the aspect ratio of the distribution. \n", "\n", "| ![Skewed data set](Skewed_pdf.png) | ![Transformed data set](Transformed_pdf.png) |\n", "|:------------------------------------------------------------:|:------------------------------------------------------------------:|\n", "| **Figure 6.1:** Skewed data set given in [Eq. 6.2](#eq:SKE). | **Figure 6.2:** Transformed data set given in [Eq. 6.3](#eq:WELL). |\n", "\n", "Let us create two Monte Carlo runs with an identical sequence of independent identically distributed variables and we use the two densities $p(x)$ and $p(y)$. If our algorithm is affine invariant, it must obey that $Y(t) = AX(t) + b$ (given initial points $X(0)$ for $p(x)$ and $Y(0) = AX(0) + b$ for $p(y)$)." ], "id": "5239a1fe0fe0fb7d" }, { "metadata": {}, "cell_type": "markdown", "source": [ "### Stretch Move\n", "Instead of using a single chain of walkers, affine-invariant algorithms use an ensemble of walkers. In this case, the updated position of each walker is based on the current position of the other walkers in the ensemble. When using an ensemble of $K$ walkers $S = \\{X_k\\}$, the position is updated by drawing a random walker from the complementary space $S_{[k]} = \\{X_j, \\forall j \\neq k\\}$\n", "\n", "Then a new position is proposed as follows:\n", "$$\n", " X_k(t) \\rightarrow Y = X_j + Z[X_k(t) - X_j]\n", "$$\n", "\n", "where $Z$ is a randomly chosen scaling factor from the distribution $g(z)$. Typically, it is proposed that:\n", "$$\n", " g(z) \\propto \\begin{cases}\n", " \\frac{1}{\\sqrt{z}} & \\text{if } z \\in [\\frac{1}{a}, a] \\\\\n", "\t0 & \\text{otherwise}.\n", "\t\t\t \\end{cases}\n", "$$\n", "\n", "This implementation of the invariant algorithm is informally called the ''stretch move''. The autocorrelation time for this method is much shorter than that of typical MCMC methods like the Metropolis-Hastings algorithm. \n" ], "id": "b75df687e6fc2d04" }, { "metadata": {}, "cell_type": "markdown", "source": [ "### The parallel stretch move\n", "Whilst using the stretch move is much more efficient than simpler algorithms like Metropolis-Hastings, the walkers are updated in series, one after another. If we could evolve them all in parallel this may significantly improve the autocorrelation time. \n", "\n", "This is not entirely feasible in EMCEE as this would violate detailed balance. Detailed balance ensures that the chain is reversible, allowing for sampling from the correct target distribution. However, an alternative is to split our ensemble of $K$ walkers $S = \\{X_k\\}$ into two subsets $S^{(0)} = \\{X_k \\forall k \\in [1, K/2]\\}$ and $S^{(1)} = \\{X_k \\forall k \\in [K/2 + 1, K]\\}$ and run these in parallel. One can observe that these two subsets are now each others complementary space $S_{[k]}$ and, therefore, we evolve the walkers in $S^{(0)}$ only by positions in $S^{(1)}$ (and vice versa) using the stretch move discussed above. Each walker in subset $S^{(0)}$ can now be evolved independent of the others in that subset. \n", "\n", "Even though the autocorrelation time does not improve significantly compared to the regular stretch move, by using this method we unlock generic parallelization, allowing for multi-threading and distributed computing, which makes this method especially applicable for modern hardware architectures." ], "id": "589c6084a8835c74" }, { "metadata": {}, "cell_type": "markdown", "source": [ "### References\n", "\n", "\n", "[Budd, T. 2024, Monte Carlo Techniques](https://hef.ru.nl/~tbudd/mct/lectures/markov_chain_monte_carlo.html)\n", "\n", "[Foreman-Mackey, D., Hogg, D. W., Lang, D., & Good-\n", "man, J. 2013, Publications of the Astronomical Society of the Pacific, 125, 306–312](https://arxiv.org/pdf/1202.3665)\n", "\n", "[Gelman, A., Carlin, J. B., Stern, H. S., et al. 2013,\n", "Bayesian Data Analysis, 3rd edn., Chapman Hall/CRC Texts in Statistical Science Series (Boca Raton, Florida: CRC)](http://www.stat.columbia.edu/~gelman/book/BDA3.pdf)\n", "\n", "[Goodman, J. & Weare, J. 2010, Communications in Applied Mathematics and Computational Science, 5, 65](https://msp.org/camcos/2010/5-1/p04.xhtml)\n", "\n", "[Nagler, T. 2021, Statistics for Astronomy and Physics Students, Creative Commons AttributionNoDerivatives 4.0 International (CC BY-ND 4.0),\n", "lecture notes, Fall 2020, Version: June 21, 2021](https://tnagler.github.io/stan-2020.pdf)\n", "\n", "[Wang, W. 2022, American Journal of Physics, 90,\n", "921–934](https://arxiv.org/pdf/2204.10145)\n" ], "id": "17b54df98cf732a8" }, { "metadata": {}, "cell_type": "markdown", "source": "", "id": "33c2ac1d5eba1542" } ], "metadata": { "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 2 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython2", "version": "2.7.6" } }, "nbformat": 4, "nbformat_minor": 5 }